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INTRODUCTION 


The Phase II Study was undertaken to investigate the closed-loop 
data reduction techniques described in the first Study Report 
and to investigate some new ideas that had evolved since that 
report was written. With the availability of the computer hardware 
at NASA/ARC, it required only the addition of the interface 
hardware between the Oscar-F Data Reader and the Computer, plus 
the associated software, to implement what was initially called 
the Phase II simulator. This simulator has now become the FILMCLIP 
System which has provided NASA/ARC with closed-loop ionogram 
processing capability since September 1968. Since that date, 
many thousands of topside ionograms have been reduced to electron 
density profiles, including, high altitude ionograms that cannot 
be processed on an open- loop basis. The FILMCLIP System will 
continue to be used in the development of new computational 
software that is described in this report. 
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ABSTRACT 


High altitude topside ionospheric sounder data is being provided 
from both the Alouette II and ISIS-1 satellites. The resulting 
ionograms are frequently so difficult that open-loop processing 
techniques are inadequate to cope with the problem. Closed- loop 
data reduction methods offer a practical solution to the problem 
of reducing ionograms from high altitude soundings, and at the 
same time provide improved capability for processing all ionograms. 

The FILMCLIP System is a Closed-Loop Ionogram Processor for the 
reduction of topside ionospheric sounder data recorded on 35mm 
Film. This system was implemented at NASA/ARC as part of a 
continuing study to develop methods of automating the data reduction 
of topside ionograms. The system has worked so well that it is 
now used nearly full time for film ionogram data reduction on a 
production basis. 

The final data reduction system, which will acquire input data 
directly from magnetic tape, is known as the TAPECLIP System, 
Software is being developed for the TAPECLIP System that will use 
much of the existing software. This software is organized 
together with new software to exploit the new closed-loop 
processing techniques that are now available. The TAPECLIP 
software is being developed and checked out under simulated 
conditions on the FILMCLIP System. 
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A number of different methods for computing electron density 

profiles from topside ionospheric sounder data are presented. 

r 1 3 1 

The method of Inverse processing is a variation of Jackson’s 1 
parabolic in log N lamination method. Lamination heights are 
preselected at equally spaced increments of true height. The 
scaled X-Trace data is curve fit so that the virtual depth, h^, 
is known for any frequency f on the X-Trace. The plasma frequency 
f N is then adjusted in an iterative procedure for each lamination 
boundary until computed points in the h' (f) plane match the curve 
fit X-Trace. 

Inverse mixed -mode processing provides a way for automatically 
combining scaled data from both the X and O-Traces to provide a 
composite electron density N (.h) profile. 

The Matrix method of ionogram data reduction uses an analytic 
model of the N(h) profile. In this method the coefficients of a 

mathematical model of the N(h) profile are adjusted to minimize 
the error between scaled and computed values of h'(f^) in a 
weighted least squares sense. 

The method of Horizontal processing provides a means for first 
order correction of the effects of horizontal gradients of 
electron density in the satellite orbital plane. The method of 
Horizontal processing recognizes that the topside sounder ionogram 
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data and the resultant computed electron density profile are 
both functions of time due to the satellite velocity with 
respect to earth based coordinates. 

New subroutines that were developed in the course of optimizing 
some of the computational software as follows : 

a. Inverse mapping to specified f 

b. Simplified solution of a matrix equation 

c. Optimal step size by Golden Section Search 

d. Matrix pseudoinverse computation 

The weighted least squares approximation method of Mallinckrodt 
was utilized in the matrix method. 

An investigation was conducted to show that values of f H , <j>, and 
h between accurately computed points could be obtained by 

S 

interpolation. The use of interpolation has greatly reduced the 
computer work load in the determination of values for these 


parameters „ 



SECTION I 


FILMCLIP SYSTEM 


1.0 INTRODUCTION 

The FILMCLIP System is a Closed-Loop lonogram Processor that is 
now in use on a production basis for the reduction of topside 
ionospheric sounder data recorded on 35 mm film. The FILMCLIP 
System was implemented in September 1968 on a minimum 
cost basis as part of a continuing study for NASA/ARC on advanced 
equipment and techniques for topside ionospheric sounder data 
reduction. Since that date many thousands of ionograms have been 
reduced to electron density profiles, including ionograms from 
high altitude low density soundings. The latter are almost 
impossible to reduce on an open loop basis due to the difficulties 
involved in identifying the portion of the X-Trace due to vertical 
propagation. 

In early 1968, the Computation Division at NASA/ARC made available 
an IBM 1800 Computer, IBM 2250 Graphic Display Unit, and IBM 360/50 
Computer to the Space Sciences Division on a second shift, 
non-interference basis, with the availability of this computing 
equipment, Astrodata proposed that a Phase II simulator be 
implemented at NASA/ARC using the existing Oscar-F data reader 
and the above computers to provide a closed-loop film ionogram 
processing capability for investigating advanced processing 
algorithms, and for reducing film ionograms on a production basis. 



The follow-on study for the development of improved topside 
ionogram data reduction methods has a number of objectives as 
follows % 

a. Improved closed-loop processing methods which take 
advantage of some of the newer optimization techniques, 

b. Relegate more of the ionogram pattern recognition task 
to computer software. 

c. Take better advantage of all the data in an ionogram 
for production ionogram data reduction. 

d. Correction of known deficiencies in existing data 
reduction methods . 

e. Use the FILMCLIP System as a simulator for development 
of TAPECLIP data reduction software. 

1.1 DESCRIPTION 

The FILMCLIP System is an on-line interactive computer processing 
system that has been operational at NASA/ARC since September 1968. 
A functional block diagram of the FILMCLIP System is shown in 
Fig. 1-1. The performance characteristics of the FILMCLIP 
System were defined by Astrodata and the details worked out in a 
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FIGURE 1-1 NASA/ARC FILMCLIP SYSTEM 
















series of conferences with NASA/ ARC and Informatics, Inc, The 
operational software for the systeki was provided by Informatics. 

A number of technical reports and papers have been written about 
the FILMCLIP System*’ 1 1 1 * to which the reader is directed for 
specific detailed information. 

The Scaling Converter in the system was built by Astrodata as a 
piece of interface hardware between the Oscar-F data reader anc| 
the IBM 1800 Computer. A functional block diagram of the Scaling 
Converter is shown in Fig. 1-2. A detailed description of this 
unit, including logic diagrams and schematics, is contained in 

f , ,i 

the Instruction Manual. 11 



FUNCTIONAL BLOCK DIAGRAM 
SCALING CONVERTER 
FIGURE 1-2 



SECTION II 


INTERPOLATION 


2 . 0 INTRODUCTION 

Significant savings in computer time can be achieved by simple 
interpolation of parameter values between accurately computed 
points, compared with direct computation of the parameters each 
time new values are needed. The feasibility of interpolating 
intermediate values of gyrofrequency f, T and satellite height h 
were investigated. The dip angle 6 changes so slowly that tests for 
interpolation accuracy for 0 were not included in this investigation 

2.1 INTERPOLATION OF f„ 

n 

In the data reduction of topside ionograms to electron density 
profiles, the X-Trace is used almost exclusively. This is 
primarily true because the X-Trace is nearly always continuous 
at the low frequency end down to the cutoff frequency, while the 
low frequency end of the O-Trace is usually missing. Since the 
reflection point of the extraordinary wave is dependent on both 
the electron density N and the flux density B of the earth's 
magnetic field, it is necessary to compute values of the earth's 

fl3] 

magnetic field each time the group path integral in eqn. (2-1) 1 
is evaluated. 
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h ' (f ) = / y’ [N (h) ,B (h) , 9 (h) , f ]dh 


where f = sounder pulse frequency 

h’ = virtual depth of reflection of frequency f 

A 

h g = true height of satellite 

h^ = true height of reflection of frequency f 
y 1 = group refractive index 

3 

N = electron density in electrons/cm 
B = earth's magnetic induction in gauss 
0 = magnetic dip angle 


It is convenient to work with the electron plasma frequency f 
and gyrofrequency f H so that eqn. (2-1) becomes 


h' (f) 
x 


= / y' (X,Y,(j>)dh 


s 


where 


X 



2 



<j> = 90° - 0 for vertical propagation 

f N = /w/raoi) 

f Tr = 2 # 8 B 

n 


( 2 - 1 ) 


N 


( 2 - 2 ) 
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Direct computation of the geomagnetic field as a subtask in the 
evaluation of the group path integral takes a significant amount 
of computer time. This computation is based on a spherical 
harmonic expansion representation of the field with coefficients 
from Daniels and Cain.^ 17 ^ 

In an effort to reduce the number of direct computations of the 
geomagnetic field, a task was set up to evaluate the possibility 
of using interpolation as a means for determining the gyro- 
frequency. The objective of the task was to provide a comparison 
between directly computed and interpolated values of f H< 

A cross section of the satellite orbital plane was set up as 
shown in Fig. 2-1 with boundaries 400 km apart vertically between 
200 km and 3000 km, and approximately 220 km apart horizontally, 
which is the average distance the Alouette II satellite travels 
between corresponding frequency markers from one ionogram to 
the next. The time between corresponding frequency markers is 
approximately 31 seconds, and the period of useful data covers 
roughly one-half the ionogram which is about 1 part in 500 of 
the orbital period. From this it is reasonable to assume that 
the velocity of the satellite is essentially constant over the 
period of the X-Trace, so that linear interpolation of f R with 
respect to time at constant height should be feasible. 
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In the vertical dimension f? H varies approximately as the inverse 
cube of the distance from the center of the earth. The objective 
here was to find at what interval it is necessary to compute f„ 

n 

so that intermediate values could be obtained with adequate 
accuracy by inverse cube interpolation. The equation for inverse 
cube interpolation between interval boundaries is 


f H 


f 


HT 


R + h T \3 

R + h J 


where 


R = ^ .f° 

r - 1 

r — V f / f 

HB 7 HT 

h = height of desired f„ 
h T = height of interval top 
h B = height of interval bottom 
h B < h < h T 

f HT ~ gy^ofrequency at h T 
f jib = gyrofrequency at hg 


(2-3) 


The detailed results of the f interpolation tests are contained 

f 3 — 5 1 

m Informatics Technical Notes . The worst 

case error for f„ in these tests was less than 1 part in 3000, 

and this occurred with inverse cube interpolation at 400 km. 

From these results we conclude that by computing the gyrofrequency 
at 400 km intervals vertically, values of fg can be obtained 
with adequate accuracy by inverse cube interpolation vertically 
and linear interpolation horizontally. 
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Where previously was computed directly for each of 3 or 4 
iterations per data point and for 20 to 30 scaled data points 
per ionogram, the direct computation for f„ can be reduced to a 

n 

maximum of 8 per ionogram for Alouette II, or a maximum of 10 per 
ionogram for ISIS-1, which has an apogee of approximately 3600 km. 

2.2 INTERPOLATION OF h 

s 

A separate part of the task, described in paragraph 2.1, was to 
check the accuracy of linear interpolation of satellite height 
with respect to time between end points computed once per 
ionogram. 


The position of the satellite as a function of time is computed 
from an orbital program in which the orbital elements for each 
satellite are periodically updated. The interpolation tests were 
run on different portions of the orbit of Alouette II including 


the case for maximum 


d 2 h 


dt‘ 


The detailed results of the h g interpolation tests are contained 

r 3 — $ i 

in Informatics Technical Notes . The worst 

case error for in these tests was less than 1 part in 2500 
in the vicinity of 530 km. From the results of these tests we 
conclude that it is only necessary to compute the position of the 
satellite via the orbital program once per ionogram. All other 
values of h can be obtained by linear interpolation with respect 


to time. 



SECTION III 


INVERSE PROCESSING 


3 . 0 INTRODUCTION 

A method of reducing topside sounder ionograms to electron density 
profiles using a lamination technique has been developed by 

-r , [ 13 1 

Jackson . For convenience this method of data reduction 
is represented in symbolic form by 

h' (f) — *" N (h) (3-1) 

where h' = virtual depth in km due to vertical propagation 
f = sounder frequency in MHz 

3 

N = electron density in electrons/cm 
h = true height of reflection in km 


Since the N(h) profile can theoretically be computed from the x and 

0-Traces, and a portion thereof from the Z-Trace, the symbolic 
representation can be expanded to 


h i (f > 

N (h) 

(3-2) 

Vlf) 

-*■ N (h) 

(3-3) 


h'(f) N(h) 

z 


(3-4) 


From the equation of the group path integral 

h 

r 

h* (f) = / p' ( X , Y , <}> ) dh 
h 


(3-5) 


Hl-1 



where y ' = group refractive index 



<(> = angle between earth's magnetic field and 
direction of vertical propagation 
h^= reflection height 


the reverse relationships can be represented by 


N(h) -^h’(f) 

(3-6) 

N(h) (f) 

(3-4) 

N(h) — h’ (f) 

(3-8) 


It follows that an N(h) profile, computed from scaled X-Trace 
data, can be cross-checked by comparing scaled and computed 
0-Trace. In symbolic form 


h^(f) -^N(h) ~^h^(f) 


(3-9) 


This is the method presently used in the FILMCLIP System described 
in Section I. 


For convenience in discussing this method, the following terminology 
has been adopted: 
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Forward Processing 


h' (f) — > N (h) 


(3-10) 


Reverse Processing N(h) — > h' (f) (3-11) 

The method of Inverse Processing, which is discussed in the 
remainder of this section, uses the Reverse Processing algorithm 
in an iterative loop. 


In the Inverse Processing method, lamination heights are 

preselected at equally spaced (50 or 100 km) increments of true 

height. N is assumed to decrease monotonically with h, and the 

N(h) distribution between laminations is assumed to be parabolic 

in log N, except the first which is linear in log N. The scaled 

X-Trace data is curve fit so that the virtual depth h’ is known 

x 

for any frequency f on the X-Trace. An iterative procedure is 
then used to find 


f . 

3 


h * ( f . ) 

* r 


h' (f .) 

_2£ 1_ 


h ' ( f . ) 
x 3 


< e < .001 


(3-12) 


When computing an N(h) profile from X-Trace data, the height of 

reflection is a function of both the electron density N and the 

gyrofrequency f H . With f N . or f,. the trial variable and h.. 

constant, the only variation in f H ^ from one iteration to the 

next is due to motion of the satellite, and that can be determined 

by simple linear interpolation with respect to time, as described 

in Section II. All of the variables on the right hand side of 

the group integral eqn, (3-5) are known, and therefore h' (f . ) 

x 3 

can be computed directly. 
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In the Forward Processing method, the gyrofrequency and 

true height h^ are both unknown variables on the right hand side 
of eqn. (3-5) . Consequently with the trial variable, a new 

value of hj must be computed in successive iterations until the 
value of f R j used in the true height calculation is the same as 

h 3] 

the actual value of f„. at altitude h . . ■ J 

H] J 


[ 1 Si 

Lockwood 1 J points up another problem with Forward Processing, 
and that is the iterative solution diverges if the slope of the 
height of reflection curve has a larger absolute value than the 
slope of the gyrofrequency curve. As far as we know now, this 
problem does not exist with the Inverse Processing method. 


Enough experimental work has been done at NASA/ ARC to demonstrate 
that the basic Inverse Processing algorithm is valid. Some of 
the techniques described in this section have not yet been verified. 


3.1 DESCRIPTION 

In Inverse Processing, lamination heights h_j are selected at equally 
spaced increments of true height; i.e. , 50 or 100 km. X-Trace data, 
scaled in the conventional manner, is curve fit so that h 1 (f ) is 

X 

defined for any frequency from f to the upper frequency limit of 
the X-Trace on the ionogram being reduced. In the TAPECLIP System, 
h ' (f) will be known for every line in the X-Trace so that curve 

X 

fitting will not be required. A simple linear interpolation will 
suffice for data points at frequencies that fall between lines. 
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For each lamination, except the first, the value of N at the 
lamination bottom at true height h^ may be approximated for the 
first iteration by extrapolation using a parabolic in log N 
equation. 


[i 3l 

The basic parabolic in log N equation , with lamination 
boundaries as shown in Fig. 3-1, is given by 


h = h j-2 + a j-l lQ 9 


N . t 

H— + b j-l lo « 
] z l 


_N 

N j-2 


a ! ' b l' b 2 = 0 


(3-13) 


h j-2' N j-2 


h . . , N . . 

D-l 3-1 


h . , N . 
3 3 


- ' 

_._V h,N 

Lamination j-1 

^ 

\ 


\ 

\ 

Lamination j 

\ 

\ 

\ . 


LOG N 

>» 


3 j~l 



-1 


Fig. 3-1. Pcirabolic In Log N Representation 

Rewriting in terms of N with N * N., the first trial value of 
electron density for each lamination is given by 

N. = N. , e z (3-14) 

J 3 J- 
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f H j + / f Hj 2 + « Nj 2 

f. = (3-15) 

3 & 

where f^ = J NN/12400 

and the computed value of virtual depth, h ’( f,), is then given 

x 3 1 

by eqn . ( 3-5) . 

In the first lamination the variation in electron density is 
assumed to be linear in log N. The trial value of N at the 
bottom of the first lamination may be estimated from the 
equation 

N 2 - N 

where a^ is derived empirically. 
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An alternative method of establishing initial values of at 
each lamination boundary was described in Section III of the 
Study Report, Part I^ 16 ^. In this method, values of are 
linearly extrapolated from the two previous ionograms in the 
same pass at constant lamination boundary levels. The 
extrapolation method should work for a large class of ionograms 
where the gradients at constant true height 




=cons tant 


(3-17) 


vary slowly throughout the pass. 


If h’(f..) < h'(f..), as in Fig. 3-2, the trial value of N., was 
X J XX J- L ]J- 

too large. For the second iteration choose 

N . 2 = N j x ( 1 + a) (3-18) 

where a “ .01 with the sign - for hT < h ' 

and + for h 1 " > h ' 


The second subscript in the above notation is the iteration count. 

If h'(f.~) > h 1 (f .„) , the trial value of N,» was too small. A 

x 32' x 32 32 

straight line thru points h 1 ( f . . ) and h ' (f . 0 ) will intersect the 

X ]I X 3 ^ 

X-Trace at f j ^ , as shown in Fig. 3-2 where the first two points 
are on opposite sides of the X-Trace, or in Fig. 3-3 where the 
first two points fall on the same side of the X-Trace. 
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FIG 


(3-19) 


For the third iteration, 



/77<f, , Hj , 


and from eqn. (3-5) a third point h^(f is computed, A segment 

of a circle passing thru the three computed values of h ' (f .) will 

x j 

intersect the X-Trace at frequency f ^ from which it should be 
possible to compute the final value of Nj from eqn, (3-19) and 
the relation 


N(h.) = 12400 f Nj 2 (3-20) 

without computing h ' (f . ) a fourth time from eqn. (3-5) . Of course 

x j 

this last conclusion needs to be verified, but if it is valid it 
will help considerably in reducing the computation time required 
to compute an N(h) profile. 


If the criterion of eqn. (3-12) is met on any of the first three 
iterations, the iterative sequence for that lamination is terminated 
and the program continues to the next lamination. 


The lamination bottom of the last lamination in the N (h) profile 
will be dependent on the scaled data point h ' ( f ) at the highest 
frequency on the X-Trace. For this case the height of reflection 
is unknown, and so the Forward Processing algorithm is required. 
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The objectives of the Inverse Processing method are as follows? 

a. Minimize the computer time required to compute an 
N (h) profile. 

B. Simplify the implementation of mixed -mode processing. 

c. Optimize the spacing between lamination boundaries of 
the N (h) profile. 

d. Make the lamination boundaries the same for all 

ionograms for convenience in interpolating and 

extrapolating values of N (h) , f„ and <j> . 

xi 
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3 . 2 STEP SEQUENCE 


The following is a description of the sequence of steps for 
computing an N (h) profile from scaled X-Trace data by the Inverse 
Processing method. The following step numbers relate to the 
processing blocks in the flow diagram of Fig. 3-4. 

1. Scale the X-Trace in conventional manner. 

2. Curve fit the X-Trace with a suitable function. 

3. Compute the plasma frequency at satellite height 

f N^ h l^ = / ~ ? 7< f l ” f Hl^ 

4. Select lamination boundaries h^ at equally spaced 
increments of true height. 

a. 50 km <_ (h^ - hj) < 100 km for step 4b. 

< 150 km for step 4c. 

b. IF (h n .LE. 1500 .) h . - m(50)km 

c. IF (h. .GT. 1500.) h. = m(100)km 

I 3 

j = 2 , 3 , . . . . , J, man integer 

5. Initialize the iteration count., ITCNT * 0 

6. Select trial value for f^j = /~W7/l2400 
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6.1 First iteration 


a. First lamination N 2 = N^e 


V h l 


where a 2 is derived empirically 

<r 

b. Subsequent laminations Nj = Nj_^e‘ 


where z 


j - /■! - 4b j-i ( Vi - V 


2b 


j-1 


a . = a . , + 2 b . , log. .. 

3 j -1 3~1 ^\N_. 


N ri 


j -2 j 


a^, b^/ b 2 = 0 


j “ 3 , 4 / • • « . t J 


6.2 Second iteration N . 9 = N (1 + a) 

where a .01 with sign - for h' < h' 

+ for h 1 > h ' 


6.3 


Compute h'ff.-.) and h'(f.~), processing Fig. 3-4 blocks 
3 -L 3 * 

7 thru 19. 


6.4 


Compute 


equation of a line thru 


h ' 

111 


h ' 

n jl 


j2 


f jl 


31 


sf 


jl 


h ' =* sf + r 

X 


the two points in step 6.3. 

slope / 

} 

intercept 
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A 


Fig. 3-4. Flow Diagram - 
Inverse Processing 





















r 


6.5 

6.6 

7. 

8 . 

9. 

10 . 

11 . 


Solve for frequency f^ at the intersection of the 
straight line and the X-Trace. 

Compute f N . - /f j3 tf j3 - f Hj ) 

/ 

/ 

Compute processing Fig. 3-/1 blocks 7 thru 19. 

x 3 l 


Increment the iteration counter 


Compute 




+ 4f 


H 


2 


ITCNT 


2 


ITCNT + 1 


Check if f . < f_ 

j max 

a. YES, the trial frequency is still on the X-Trace. 
Compute h'(f-) on branch to Fig. 3-4, block 10. 

X J 

b. NO, there is no more X-Trace data. Branch to 
block 22, Fia . 3-4. 


Compute Tj == t(fj) from time and frequency marker 
calibration data. 


At constant h 


with i = 1 , 2 , . . . . , j is computed 


‘i ' Hi 

by linear interpolation with respect to time, as 
described in Section II, to compensate for the change 
in satellite position. 
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12 . 


13. 

14. 

15. 

16. 


In the first lamination the slope of the N(h) curve, 
which is linear in log N, is given by: 


h. 


a. 


log 


N, 

N, 


Eqn. (2) 


[ 18 ] 


where h^ = satellite height 


As the computation of N(h^) progresses for each 
lamination boundary, the change in satellite height 
can be determined by linear interpolation with respect 

to time, given the slope of the orbit for that ionogram. 

Ah. th 

n s. The new value of N (h. T.) for the j lamination 

At" 1 3 


is computed from 


h lj" h 2 


N lj ‘ N 2 e 


Compute 


'Hi 


f . 
3 


2 , 


n 3 


Compute 


Y. 

i 



i = 3, 4, ] 

a^ , , b£ = 0 

Eqn. (3) , (4) , [l8] 
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17 . Compute h'(f.) using eqn. (7), 

X J 


[ 18 ] 


18. Compute error ratio and compare with limiting value 


!h x (f 1> - 




h’ (f .) 
x ' y 


< .001 


19. Compare iteration count with integer 3. 

a. If <, continue next iteration, step 6. 

b. If compute N(hj), steps 20 and 21. 

20. Pass a smooth curve (segment of a circle) thru points 

Nctfjk)' k = 2 ' 3, intersecting the X-Trace at fj^. 

21 . Compute f Hj - f”“) 

N (h j) = 12400 f Nj 2 

22. Compute N(hj) by Forward Processing, where hj is the 
bottom of the last lamination. Output N(h-), 

i = 1, 2, j 

j « J 
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3 . 3 MIXED MODE 


The useful portions of the O and Z-Traces can be used in the 
Inverse Processing method as well as X-Trace data as described in 
in paragraph 3.1. 


X, Of, and Z-Trace data are scaled in the conventional manner and 
curve fit so that h^(f), h^ (f ) , and h^ (f ) are defined for all 
frequencies on the useful portions of their respective traces. 
Mixed-Mode Inverse Processing begins with X-Trace data as in 
paragraph 3.1. At each lamination boundary at the point in the 
program where the sounder transmitter frequency f ^ is computed, 
as in eqn. (3-15), the corresponding frequencies for the 0 and 
Z-Traces are also computed. 


'*3 


f.+ /f . 2 + 4f ? 
j / Hj x Nj 


f . « f„. 
03 N3 


f . = f . - f„. 
Z3 xj H3 


(3-21) 

(3-22) 

(3-23) 


A comparison is made with the curve fit 0 and Z-Trace segments 
to see if there is an h^(f) at frequency f . or an h^ (f ) at 
frequency If . . For this example, assume that an h^(f Q j) does exist. 
The Reverse computation 


N (h 


j> --^v 


(3-24) 
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gives the first point for lamination j on the h' (f) plane. The 
iterative sequence continues as in paragraph 3.1, except for the 
third iteration where 




(3-25) 


from which a third point h^f^) is computed. A segment of a 
circle passing thru the three computed values of h^(fj) will 
intersect the 0-Trace at frequency f^, from which it should be 
possible to compute the final value of from eqn. (3-25) and 
eqn. (3-20) . 


To compute IQ (h^) with respect to curve fit Z-Trace data use eqn. 
(3-23) and the Reverse computation 

N(h.) — * h' (f .) (3-26) 

J ^ J 

At any lamination boundary h^, it would be possible to compute an 
N(hj) with respect to each of the X, 0, or Z-Traces , assuming the 
traces have been scaled at the related sounder pulse frequencies. 
The value of N(h^) would be slightly different for each trace 
and would depend of course on the delay contribution of the 
previous laminations. 


On the basis that the O-Trace is usually thinner than either the 
X or Z-Trace, the program logic for Mixed Mode Inverse Processing 
might be to compute N(h^) with respect to the 0-Trace whenever 
possible, otherwise use X-Trace data. When processing with 
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respect to the 0-Trace, corresponding h’ (f .) 

x j 

computed to assist in defining the vertical 
of a difficult X-Trace. 


points could be 
reflection portion 
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SECTION IV 


MATRIX PROCESSING 


4.0 INTRODUCTION 

One of the problems an operator faces when processing ionograms 
with the FILMCLIP System occurs after one processing iteration 
represented by 

h ;(fj) -*N(h.) h; (fj) 

If the computed O-Trace data points h^(fj) don't match closely 
enough the identifiable O-Trace on the ionogram, the operator 
must change some of the scaled data points on the X-Trace in an 
effort to modify the N(h) profile so that the computed h^(fj) 
data points will match the O-Trace more closely. This is really 
a very complex, multi-dimensional, non-linear problem that a 
human being is not equipped to handle, except by trial and error. 
With the FILMCLIP System an operator can make adjustments in the 
positions of the scaled data points on the X-Trace and observe 
the results referred to the O-Trace in successive iterations 
until a suitable match is achieved. With experience, operators 
become quite skilled in making adjustments in this iterative 
loop, but it still remains a time consuming and inefficient 
procedure . 
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The matrix method provides a logical means whereby the computer 
can use all the normally available information in an ionogram 
to compute an N(h) profile in a weighted least squares sense* 

The computation time to compute an N (h) profile by the matrix 
method will increase substantially over the time required for 
either the Forward or Inverse methods. However, it will take 
the operator out of the iterative loop for closed-loop 
processing so that there may be an overall reduction in elapsed 
time. In addition, the matrix method offers some unique advantages 
not available in other methods. 

The matrix method of ionogram data reduction uses an analytic 
model of the N(h) profile rather than the lamination model of 

f 1 3l 

Jackson . The concept behind matrix processing is to set up 
a mathematical model of the N(h) profile and then adjust the 
coefficients of the mathematical model to minimize the error 
between scaled and computed values of h ! (f^) in a weighted least 
squares sense. 

A set of scaled h ’ ( f data points from an ionogram is the input 
prescription vector to the matrix program. The data points can 
be scaled from both the X and 0-Traces and even the Z -Trace when 
it is available. A weighting coefficient is assigned to each 
data point using an arbitrary scale from 1 to 10. Good data 
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points are weighted higher than questionable data points so that 
the resultant N (h) profile is influenced more by the data points 
in which the operator has greater confidence and to a lesser 
extent by data points with a lower weighting factor. 

The particular model used for matrix processing is a ratio of 
polynomials of the type used in synthesizing electrical filter 
transfer functions. It is perhaps one of those fortuitous 
circumstances, but when a set of N(h) values were input to a 
slightly modified Astrodata proprietary program for filter 
synthesis, a curve fit within 1% was achieved the first time the 
right number of poles and zeros were specified. 

With the matrix method it is desirable to minimize the number of 

coefficients required to model the N(h) profile, because some 

3 

matrix operations on the computer go up as m where m is the order 
of a square matrix. The roots of the numerator and denominator can 
be plotted as poles and zeros in the complex s-plane. A plot of 
the poles and zeros for one model is shown in Fig. 4-1. In this 
model there are four conjugate complex poles, two real poles, and 
one zero at the origin. Other models investigated have used five 
or six complex poles with from two to five zeros at the origin and 
no real poles . The number of adjustable coefficients required for 
these models varies from 10 to 12 . Kinkel ^ 2 1 ^ proposed a model 
with 5 variables which turned out to be an insufficient number . 
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The choice of the particular zero/pole representation for an 
N (h) profile depends on a number of factors, such as the ratio 

of, minimum to maximum electron density and the shape of the curve. 

f 

A (Little experience in curve-fitting a variety of N (h) profiles, 
using Program NHMODEL in Section 4.4, will provide the, basis 
for selecting the initial values for coefficients of the model 
used in the matrix method. 


I 

I 

t 

} 

I 
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4 . 1 DESCRIPTION 


The general form of the analytic model used in matrix processing 
is shown in eqn. (4-1) . 


N(s) 


ps 11 

(s+o^Ks+o ) 2 + B m 2 ] 


(4-1) 


where s 
h 
n 

P 

a 

e 

j 

z 

m 

n 


j (n-h) 

3 

true height above mean sea level m km x 10 
a constant > maximum satellite height in same 
units as h 
scale factor 

real component of pole position 
imaginary component of pole position 

/=i 

0 , 1 , 2 index of real poles 

1, 2, . M index of complex poles 

1, 2, . 5 index of zeros at origin 


Eqn. (4-1) is a representation of an N(h) profile in the complex 
s-plane using two real poles, four conjugate complex poles, and 
a zero at the origin. A representative distribution of poles and 

j 

the zero are shown in Fig. 4-1./ In this representation, only 
the magnitude of the function is considered. The denominator of 
the function for a given value of s^ on the positive imaginary 


I 

i 
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i 







axis, is computed as the product set of all the vector magnitudes, 
| p k l ^ . The numerator of the function is the vector magnitude z 
along the imaginary axis from s^ to the origin multiplied by the 
scale factor p. 

Then the value of points N(s^) is given by eqn. (4-2). 

N (s i ) 


P z. 


n 


K 

TT 

k=l 


'k i 


(4-2) 


where s^ = j (p-h^) as in eqn„ (4-1) 

1 *™ 1, 2 , Baas, X 


Define a column vector 


a = 



a 

a 

a 

a 


1 

2 

3 

4 


a 

a 

a 

0 


1 

2 

3 

3 



a K-l “ °M 
a K = 0 M 


(4-3) 


Ml of the variables a and 3 in eqn. (4-1) are represented by a. 
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In this notation 


I = number of scaled data points in the 
prescription vector 

K = number of degrees of freedom in the 

mathematical model of the N(h) profile 

An initial set of values of the a vector is established to 
represent a trial function of N (h^ , a) for a set of scaled h^' (f^) 
data points from an ionogram, where i = 1, 2, . ..., I. The 
hf' (f^) will be a composite set of points scaled from both the 
X-Trace and 0-Trace, and even the Z-Trace if it is available. 

We wish to approximate the set of scaled data points by a 
vector H(f i , a) in the h' (f) plane which is a mapping of points 
from the approximation N(h^, a) in the N(h) plane as shown in 

f 2 3 I 

Fxg. 4-2. These points are computed using Reverse Processing 

as defined in Section III. The value of h^ for each point 
N ( h a) is computed by successive approximation in an iterative 
loop so that 



This iterative loop is described in paragraph 4.2. 
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> 




Fig. 4-2. Mapping of Corresponding Points in 
the h' (f) and N(h) Planes 


The differential of the function H(f^, a) can be approximated 
by a small increment on each a v as in eqn. (4-5) . 

3H(f. , a) 

Ao k (4-5) 

i 1/ 2 , .«../ X 


K 

AH(f., a) = X) 


Eqn. (4-3) expressed in gradient notation becomes 


AH(f i , a) = «)% 


(4-6) 


1^ 2 1 . . » • f X 


where 


V H(f . , a) 
a l — 


9H(f i , a) 
3al 


3H(f i , a) 


3 a 


K 


Aa « ( Aa -^ , ...., Aa^) 


T 
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We wish to find 


AH(f i/ a) 5 h'(f i ) - R(f ± , a) - AHff^ a) = 0 


1 3 1 , 2 y 


I 


(4-7) 


Rearranging terms and normalizing by dividing both sides by 
H(f i# a) 

h' (f . ) - H(f . , a) AH(f . , a) 

i__ i — Z_ . _i — — (4-8) 

H(f i , a) H(f i# aj 


Combining eqn. (4-6) and eqn. (4-8) and expanding in matrix form 
gives 


a ll a 12 a lK 

a 21 a 22 a 2K 


a Il a I2 ‘ ’ a lK 


Aa, 


Aa. 


Aa 


K 


where e 


h« (f ± ) - H ( f t , a) 
H(? a , a) 


a ik 


9H(f^, a) 

H (7^, a) 9a k 


i = 2, 3 , 

k = 1, 2, 


I 

K 


(4-9) 
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h' (f x ) = 0 


For i = 1, 


where 


f l = f KS 


consequently e^ and k = 1, 2, . » . . , K, must be redefined in 
terms of frequency because AH has no meaning at satellite height, 
For this case 


e l " 


f x (h x ) - F(h r a) 
F(h x/ a) 


lk 


3F(h 1# , a) 
F(h 1 , a)3a k 


In compact matrix notation 

E = A Aa 


(4-10) 


In eqn. (4-9) there are I equations with K unknowns. With 
I > K there are more equations than unknowns since typically 

20 £ 1 <. 30 
and 10 < K < 12 


A least squares fit solution can be obtained by first 
multiplying both sides of eqn. (4-10) by A^, which gives 

A T E - [A T A]Aa (4-11) 

where T = the transpose 
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Then a least squares estimate of the adjustment vector Aa is 
given by 

A a « [A T A] A T E (4-12) 

We wish to assign a weighting coefficient, w^^ , to each scaled 
h' (f^) data point in order that good data points exert a greater 
influence on the least squares approximation of N(h, a) than 
questionable data points. Define a diagonal matrix of weights 
as follows s 


w 


11 


w 


W = 


22 0 
W 33 


w 


II 


(4-13) 


The weighting matrix W is incorporated in the least squares fit 
solution of Aa to give a weighted least squares adjustment vector 
as follows: 

Aa « | A T Wa] A T WE (4-14) 
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For each iteration then 


where 


a 


n +1 


a 


n 


/N 

+ Aa 



n = iteration count 


(4-15) 


which gives a new set of coefficients for the analytic model 
N(h, a). The iterative 9 yd e continues until 


A a, 


a. 


< e < .001 


max 


(4-16) 


A more direct measure of how well the computed data points 
H(f^, a) match the scaled data points h'(f^) might be obtained 
by defining a cost functional J as the sum set of the weighted, 
normalized error squared. 

I 

J = X) d i (4-17) 

i=l 

where D = [d^, d^r ••••» d.j.] T = 


The iterative cycle would then continue until 


T n +1 T n 

J ~ U 


,n 


< e 


where 


n = iteration count 


(4-18) 


A flow diagram showing the major processing steps in a program 
for the matrix method is shoton in Fig. 4-3. 
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4 . 2 STEP SEQUENCE 


The following is a description of the sequence of steps for 
computing an N(h) profile from scaled X-Trace data by the Matrix 
Processing method. The following step numbers relate to the 
processing blocks in the flow diagram of Fig. 4-3. 


1. Scale the X, 0, and Z-Traces of an ionogram in 
conventional manner. 

2. Assign a weighting coefficient to each scaled 
data point. 

3. Establish initial values of coefficients for a. 
Program NHMODEL has been developed for doing this 
and is described in paragraph 4.4. 

4. Compute values of h^ corresponding to frequencies f^ 
of scaled data points h'(f^) such that 



r 2 3 I 

This is accomplished by a one-dimensional search 1 1 
using successive approximations and is described in 
paragraph 4.3. This makes possible comparison of 
computed and scaled data points at the same frequency 
in the h '( f) -plane. 
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Fig. 4-3 
FLOW DIAGRAM - 
N (h) Processing 


Matrix Method 














5. Compute data points H(f^, a) in the h ! (f) -plane using 
the Reverse Processing algorithm of Section III. 


N(h ± , a) — > H ( f ± , 




a) 




h i (f > 

h;<f) 

5JT?) 


The electron density profile N(h^, a) for this computation 
is obtained from the mathematical model of eqn. (4-1) . 

The computation of N(h^, a) is handled by SUBROUTINE CALC 
and is described in paragraph 4.4. 


6. Compute the elements of discrepancy vector E as 
defined in eqn. (4-7) . 


7. Compute the elements of matrix A. This is the matrix 
of partial derivatives of H(f^, a) with respect to a, 
with each element a^ as defined in eqn. (4-7) . The 
term 3 a^. is approximated by perturbing the elements of 
a by 0.1%, one at a time. 

8. Compute the elements A a. of the adjustment vector Aa, 
•which is the solution of matrix eqn. (4-12) . This is 

handled by SUBROUTINE MXV which is described in 
paragraph 4.4. 
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9. This step is a decision block to determine whether or 

not the latest adjustment vector will make a significant 
change in the coefficients of the mathematical model 
of the N(h) profile. If the maximum absolute value of 
all the change ratios of vector a are less than e, 
written as 

% 

Aa, 

< e, k = 1, 2, . , K 

^ max 

a. NO, branch to Step 10. Repeat steps 4 thru 9. 

f 

b. YES, a weighted least squares solution for N(h^, a) 
has been obtained. 

A 

10. Add the adjustment vect6r Aa to the state variables a and 
repeat steps 4 thru 9. 

11. Output N(hj, a) at equally spaced increments of true 
height h ^ . 
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4,3 COMPUTATION OF h ± 

The matrix method of computing an N(h) profile requires that 
comparisons of h'(f) and H(f^, a) be made at the same frequency. 
Since only portions of an 0-Trace are usually available, and we 
would like to avoid the extra step of curve fitting a scaled 
X-Trace, a method is described in this section whereby values 

ft* 

of true height h^ are computed such that computed data points 
H ( f ^ , a) will occur at the same frequencies as the corresponding 
scaled data points h' (f^) . 


The frequency, f^, of a computed data point H(f^, a) for the 

X, 0, and Z-Traces is a function of plasma frequency f^, and 

for the X and Z-Traces is also dependent on the gyro frequency 
f 1 ;1 1 

as follows L J : 


* , / f 2 

"H / H 


(4-19) 



where N = ionospheric electron density 

3 

(electrons/cm ) 

B « induction (gauss) earth magnetic field 
f = frequency in MHz 
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Since N and B are monotone decreasing in h, and f^ are also 
monotone decreasing, and therefore the computed frequency f^, 

computed from eqn. (4—19) through eqn. (^4-21), is monotone decreasing 
in h. This property makes it possible to use a technique of 
successive approximation for finding h^ such that 



< e < .#01 


(4-24) 


For the general case 

h. < h = satellite height (4-25) 

i — s 

so that the first trial value of h^ will be 

h il - "T (4 - 26) 

If f^ > f^, then h ^2 > On the first iteration, in this 

example, h^ is in the upper half of the interval h g to h = 0. 

On successive iterations hu is bounded by 1/4, 1/8, etc., of the 
total height interval until the conditions of eqn. (4-22) are met. 

At each trial value of h.^, f N is computed from eqn. (4-2) and eqn. 
(4-22) . 

„ n 

z . 

N(n-h.) « P (4-2) 

IT l p kii 

k=l 
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where the variables are the same as previously defined. The 
value of f fl at h^ is obtained by inverse cube interpolation 
(refer to Section II) at the satellite position, which is known 
from t(f^) and the topside sounder orbital parameters. 

The following is a brief description of the various processing 
steps in sequence as th§y relate to Computing h^ prior to obtaining 
H (f^, a) by Reverse Processing. The following step numbers relate 
to the processing blocks in the flow diagram of Fig. 4-4. 

1. Input f . , t(f.), h , and trace designator (X, 0, or 

X X S 

Z-Trace) 

2. Compute trial value of h^ 

3. Compute f N , f H , f^, Ah. The value of f H (h) is obtained 
by inverse cube interpolation vertically with respect to 

t 

true height and linear interpolation horizontally with 
respect to time. 

f, - 7 . 

4. If ~=~j — - < r. 

1 i 

NO, go to step 5. 

YES, h^ = h. RETURN to calling program. 
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h = h + Ah 


7 



Fig, 4-4 
Flow Diagram - 
To Compute 

for Matrix Processing 








5. If (f i - f i ) < 0 
NO, go to step 6. 

YES , go to step 7 . 

6. Compute next trial value, h ® 

7. Compute next trial value, h 


h - Ah 


h + Ah 
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4.4 PROGRAM NHMODEL 


The Matrix method of reducing topside sounder ionograms to 
electron density profiles requires an initial estimate of N(h) 
as represented by a mathematical model. In other words, an 
initial set of values for a are required. This is Step 3 in 
the flow diagram of Fig, 4-3. 

While the program for the Matrix method is under development, 
it would be well to start with an initial function N(h, a) that 
is reasonably close to an actual N(h) profile computed by 
either Forward or Inverse Processing. Program NHMODEL was 
developed in order to curve fit a variety of different N(h) 
profiles during the investigation of ways to adequately model 
an N (h) profile. Some of the subroutines such as CALC and MXV 
can be used directly in the overall Matrix program. Other 
subroutines such as COST and STEP would be modified for use 
in the flow diagram of Fig. 4-3. 

Program NHMODEL was developed using a time share computer 
terminal. The program is written mostly in FORTRAN IV with 
some variations that were permitted and others required by 
Tymshare FORTRAN IV. Differences in the variable names in the 
program and the nomenclature in paragraph 4.1 arise due to the 
fact that Program NHMODEL is a revision of a standard Astrodata 
program and it was convenient to minimize changes in the 
variable name structure. 
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A complete listing of Program NHMODEL is included starting on 
page IV-34. The numbers in the left hand column of each page 
are line numbers that are used in the time share mode of 
operation. They are convenient for reference but are not part 
of the actual program. 

A flow diagram of program NHMODEL is included starting on 
page IV- 33a. The results obtained from curve fitting the N(h) 
profiles from four different ionograms are shown in Tables 4-1 
through 4-5. 

Input data is set up in the DATA statements, lines 138 through 144. 

NCP = number of complex poles 

NjRP = number of real poles 

NRZ = number of zeros at the origin 

The last two elements of array ZK are the scale factor p, and the 
change of variable constant n . Scale factor p is modified on 
each iteration at lines 580 and 824, but n remains constant 
throughout the program. In line 162, K is defined as the number 
of degrees of freedom in the N(h^, a) model. 

Array IIN is a two-dimensional array which contains the input 
prescription vector, which is the N (h^) profile to be curve fit. 

The (HN(I,1), I = 1 , KL) are values of true height in units of 
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MNEMONICS 


NCP 

NRP 

NRZ 

ETA 

ITCNT 

CK 

EMAX 

QS 

J 

PZK ( J) 
ZK(J) 
I 

HN (I, 1) 
HN (1,2) 
AMP ( I ) 
E(I) 


= Number of Complex Poles 
= Number of Real Poles 
= Number of Zeros at Origin 
= Constant for Change of Variable s = n -h 
= Iteration Count 
= e in eqn . (4-33) 

= Cost Functional, |E(I)|max 
= Step Size in eqn. (4-42) 

= Index for ZK Array 

= Initial Pole Positions of Analytic Model 

= Final Pole Positions of Analytic Model 

= Index of Input Prescription Vector HN 

3 

= True Height m 10 km 

-4 

= N (Electron Density x 10 ) 

= Value of N from Analytic Model 
= Curve Fit Error, £n[HN(I f 2)/AMP(I) ] 


Table 4-1. Mnemonics for Tables 4-2 through 4-5 



) • 




\ 

s: 

x 

\ 



in 

sr 

*0 


Si 



KJ 

to 


sx 

o — 

— 

S: 


RUN 

CO 

CO 

o 

II 

o 

s 

Ui 
>< K5 
<L 

s: to 

r^ 


UO P- 

a 

X 

_* 

O 

tr 

m 


sz 

in 

1 


« 



in 


u 


\ 



m 



<r 

ii 

H 

h~ 2 


32 

C X 

a. 

o 

tv;; ’ • c 

CJ 

f- 

/ 

t 

> 

s: 

f**-H 


^t^incu^^co — op-oc 
^-'fOsr<\]*n<£ f ODOCV3 0%'0 
Sscect^r- oin-inm^fo 
N^osn WOCGOO OOJ 
o co co o cs in m co co 

k> co — > *o — * <r 


o o 
n- o 

0> O 

r- o 
m o 

9 © 

<3- -q- 
<7 
<T 
00 
02 


tsocoooooocooo 

^oooooooooooo 

5COOCOOOCOOOOO 
N^VOOOOOOCCCO o 
do^'^c*. o o mm w^oo 

««<&990®99«*9 

t-C * CO — K> ^ — * •q' 


Co to c ir Vt r- o v o — Co 



CO 

<r ^ 

« 

"5T 


o oo m 

*3* O "™* 0, ■*— * fO CO 



0 

VOfcOCOf^ — 

r- r- cc — , t£> 

00 

— K) cj <«e &, — 

«— * 


#— t K5 

h^CJOO 



•<roP~ecop-ccocJ 

N-/ 

o 

o © o o o 

o o o o o 

o o c o c c o 

o o 

UJ 

o 

o o 

o o o 

o o o o o 

o o o o o o c 

o o 


<® 

9 « 

• 

» • 

9 

• • m ® 

9 

• • ® • ® ® 


* 



I 1 

8 

1 8 




1 i i 

8 

« 

/N 

C c 

cj cc -sr 

c?o m 


o r- 10 0 , 

mo; vcco^o 

sr 


X) 

>o — 

t» Cs3 O 

ho fo co <r m 


o co — ■ . cj r» c%i cj r- 


o 

in o in oj o 

m © to oc co 

r- 

iTvctxjc-r'ti-^K: 

X 

ft 

o co m in c- 

00 

— so in — 

ocr- oo^r^ k)oc 

m 


3o 

SO FO 

tO »o J-O 

FO 

<r <r <r in 

in 09 n** a . o co in 

c 


«£ 

» 

. -» e 

<9 

« « 

9 

® © © * 

9 

® 9 ® ® ® 9 


* 









■""* 

CO 

JO 


a 

o o 

o 

o c 

o 

o o o © 

O 

o o o o o e 

c 

o 

CO 

o 

o o a 

o o 

o 

c o o o 

c 

o o o o o c 

a 

c 


to 

m oc 

m co o 

o 

— «a* cv c 

o 

^ co o a cj o 

w 

r- 

1 — 1 

V ■ 

o — 

hO 

m r- 

c. 

— . *o m — 

& 

r- o, in co o 

c 

<r 


CO 

fO fO to 

K> IO 

K3 

-q- <r tn 

in 

C?C CO 

o 

— j 

52 

0 


9 

« 9 

9 

0 9 9 9 

9 

® © * * 9 * 

® 

» 

X 








^4 I «— 

Co 

tO 

/-V 

o 

o o 

o 

o o 

o 

o o o o 

o 

o c o o e c 

o 

c 

— * 

o 

o o 

CO 

o o 

o 

o o o c 

o 

o c c c c o 

c 

o 

«* 

CO 

o o 

o 

o o 

a 

o o o c 

o 

o e o cc c 

c 

cv 

F— ‘ 

o 

a c 

o 

o o 

o 

o c o o 

o 

c c c c- o o 

o 

o 


0 

a n 


it , q* 

to CO — o oc 

09 

<r cj o c o r- 

OS 

in 


9 

• 9 

• 

• • 

• 

• • « « 

9 

9 O « ® # * 

9 

© 

X 

CO 

CJ CJ 

CO 

CO CM 

CO CO CO Co — 


“ *“~ 4 



»-*-* 

_ 

Cv] rr 


m ov 

r- 

0 o o — 

Co 

kd sr in o. n 0 

c 

c 


- cv 


u: a 
cn 

X c 
<T CJ 

0. F, 


Table 4-2 » /HN1/ 65-340-003823 Alouette II Pass 81 at Ottawa 



1-13-70 


RUN 2 


1 400 


/HN4/ 


/DATA.5/ 

t 

2 4n >RUN 
SC° - 6 MRP 


0 NRZ 


3 ETA = 5 


TCNT 

CK 

EMAX 

GS 

1 

1 

,4405434 1 


o 

] 

6. PI 630 IE-02 1 

,0669977 

3 

.3 

4 .6*8824€E“02 

.6331 4002 

A 

<>E-Q2 

2.873801 3E-02 

1 .3239071 

e 

2.7E-02 2 * 395 4022E' 

-02 .69599374 

6 

IE-02 

2 .1 3*7658E-02 

.380101 41 

7 

1 E-02 

1 ,qo2«5c>24E-02. 

.21983965 

rr 

1 E-02 

1 .OP35«54E-02 

,1 1 457904 

C 

1 E-n? 

1 .O808968E-02 

7.2008353E-02 

f 0 

1 E-02 

! .O7QS544E-02 

4.5704904E-02 

! ! 

1 E-02 

1 ,Q7cf«s.oo2E-02 

2 . 3237761 E-02 


J 

P7K (J) 

ZK (J) 

1 

.06000 

.05240 

2 

7 .00000 

6.58393 

3 

.0*000 

.0641 8 

A 

6.1 1000 

5 .50347 

5 

.05000 

.05270 

f 

4 .€<*000 

4 .65728 

7 

.30000 

. 42664 

>" 

4 .7 "000 

4,7 i € A •* 

o 

.50000 

i .32791 

10 

3 .2.QOOO 

2 .9*960 

1 1 

1 .02000 

.68869 

12 

1 .75000 

1 .76945 

1 3 

1 .00000 

124084.43000 

1 4 

5 .00000 

5,00000 


I 

HN( I , 1 ) 

HN(I,2> 

A MP ( I ) 

E(I) 

1 

2.00066 

.4121 1 

,40895 

,00771 

o 

2.04012 

.43946 

.43731 

,00491 

3 

1 .0 3° 3o 

.50675 

.50802 

-.00251 

4 

1 .73615 

,67384 

.67809 

-.00628 

6 

1 .53350 

.8371 1 

,93364 

.00371 

6 

1 .352*5 

1 .28677 

1 .29261 

- .00453 

7 

1 .167 40 

1 ,00108 

1 .01385 

- .00669 

p 

1 .04240 

2.64*45 

2 .61005 

,01 46.4 

0 

,<*4030 

4 ,7f *4j 

4 *68601 

.00680 

10 

.65*50 

0 ,Q6' 7 43 

SO,! 4657 

-.01 78t 

1 1 

.65923 

17.02012 

1 8.10745 

-.01040 

! ° 

.4**1 * 

88. 1 0641 

2.8.28757 

-.00642 

1 3 

.44505 

41 .50047 

40 ,92579 

,01394 

1 4 

.41461 

56.1 8123 

55.56875 

*0 1 096 

1 5 

.3*580 

74.68321 

76,17590 

-.01979 

16 

.36753 

02.61606 

91 ,54081 

.01168 


STOP 
(3380 ) > 


Table 4-3. /HN4/ 68-122-165346 Ionogram Recorded at Ottawa 
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-70 RUN 2 1120 /HN5/ 
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/ P A T A 5 / 1-13-70 RUN 5 1550 /HN6/ 

t 

>LI ST 13R;139 

1 3 a DATA NCP, NRP, NRZ /6, 0, 3/ 

1 3 C DATA (ZKCK) , K = 1,14) / .06 , 6.6 , .0 7 , 5 . 5 , 

.056,4.7, .45,4*75,1 .3 , 3 .2 , .85 , 1 .7 , 1 . , 5 . / 

>RUN 

NCP r 6 NRP = 0 NRZ = 3 ETA = 5 


NT 

CK EM AX 

1 .1Q462359 1 

OS 



1 2.6202655E-02 

! ,0896104 


.3 2,0441 304E-02 

.8997981 6 


° E-02 2.0446496E 

-02 -3. 

0656876E-0 

J 

P7K ( J) 

ZK ( J) 


1 

.06000 

,05984 


o 

6.60000 

6,62509 


3 

.07000 

,06633 


4 

5.50000 

5 ,50698 


5 

,05600 

.05992 


6 

4,70000 

4 ,68798 


7 

.45000 

.34404 


o 

4.75000 

4,829 79 


o 

1 .30000 

1 .42056 


! 9 

3.20000 

3.27129 


1 1 

.*’5000 

.,9 \ 5 50 


12 

1 .70000 

1 ,64677 


1 3 

1.00000 160412.94900 


! 4 

5 .00000 

5.00000 


I 

HNCI.l) HN( 1 ,2) 

AMP (I) 

E(I) 

1 

2.160^4 ,33260 

.33223 

.00112 

g 

2.06682 . 38499 

.38419 

,00208 

3 

1,86519 .52641 

.52962 

- .00608 

4 

1.61567 ,20906 

,80973 

00083 

5 

1.32034 1.25600 

1 .24366 

.00987 

6 

1.10733 1 .22959 

1 .88959 

- .00000 

7 

1.06615 2.63100 

2,62495 

,00230 

P 

.02630 3.93651 

3 .93452 

,00051 

c 

,7**0? 6.22601 

6.3491 5 

-.00999 

10 

.66704 10.62651 

1 0 .58890 

,00355 

1 ! 

,56003 18,88184 

1 8.621 87 

.01386 

12 

,40175 88 ,°6? 5? 

29 .12638 

-.00564 

I 3 

.44155 42.40472 

43 .28067 

- .020 45 

1 A 

,40927 57,25030 

58,1 5866 

-.01504 

1 * 

,3050° 75.7556R 

74.22366 

.02043 

1 6 

. 3 6 1 2 o 96,451 32 

95 .30206 

.01199 

1 7 

.33612 11 R. 03530 

120 .87486 

-.01618 

! ° 

,32266 1 31 ,76640 

131 .98951 

-.00169 

! *5 

,3001 1 1 39.1 31 70 

137,71923 

,01020 


ST^° 
(•®3<>0 )> 


Table 4-5. /HN6/ 68-122-165450 Ionogram Recorded at Ottawa 
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3 

10 km. The HN(I,2) are the related values of electron 

4 

density divided by 10 . KL is the number of elements in the 
prescription vector. 

The program is initialized in lines 148 through 172. At lines 174 
and 176, the ZK array is stored in PZK for comparison on final 
printout. 

At lines 178 and 180, the change of variable computation takes 
place . 

s. - g - h. from eqn. (4-2) 

l 1, 2 , ...., 1 

Subroutine CALC, which is first called at line 188, computes 
N(s^), eqn. (4-2), as AMP(I), using the values of ex in array ZK. 

In line 580, the scale factor p is computed initially so that 
in line 584 

AMP (10) = HN (10,2) (4-27) 

This forces the two curves N(h) and N(h^, a) to cross at the 
tenth lamination boundary. 

In Subroutine COST, which is first called at line 196, the 
discrepancy vector E(I) is computed in line 816. This is similar 
to vector E in eqn. (4-10) . In the same DO loop the sum set of 
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E(I) is accumulated so that the mean value can be computed in 
line 820. The mean value of E(I) is removed in line 830 which 
is equivalent to changing the scale factor P by a term called 
EXM in line 822. The scale factor correction takes place in 
line 824. On the last iteration, the scale factor correction 
is applied to AMP (I) in line 314 for printout. In the DO loop 
lines 828 through 836, the maximum absolute value of E(I) is defined 
as EMAX . 

In Program NHMODEL, EMAX is used as the cost functional and the 
program is structured to minimize EMAX, although it minimizes all 
other values of E(I) as well. Beginning at line 838, if EMAX > .5 
a different scale factor is computed as ER. This is used to 
scale down the discrepancy vector E(I) in line 844 so that the 
maximum value of E in the matrix equation does not exceed 0.5. 

This prevents a divergent condition from developing if the 
initial values of ZK are not very good. 

Statement 100 at line 198 is the reentry point for the major' 
iterative loop and is called at line 298. There are three exit 

conditions from this major loop at lines 202, 208, and 209. 

/ 

The elements of the matrix of partial derivatives D(I,L) are 
computed in Subroutine GEND, which is called at line 216. 

These are similar to the elements a^ of eqn. (4-9), but the 
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D(I,L) can be computed directly because N(h, a) is an analytic 
function of a and is the model of the curve to be fitted. The 
elements are computed as the partial derivative of the natural 
logarithm of the function N(s, a) with respect to a as given by 


D(I,L) = 


9 An N(s.) 
3 a, 


(4-28) 


where the notation D(I,L) is similar to a^ in eqn. (4-9) 


So far in the program we have developed a matrix equation 

E = D(DK) (4-29) 

where DK is the program name for the adjustment vector Aa- The 

similarity with eqn. (4-10) is obvious. Both sides of eqn. 

T 

(4-29) are multiplied by D in lines 224 through 240. In the program 

B = D T E (4-30) 

and A = D T D (4-31) 


Since A is a symmetric matrix, only the terms on and above the 
main diagonal are computed directly. Terms below the main 
diagonal are filled in at line 254. 


Under some conditions matrix A can be poorly conditioned so that it 
cannot be properly inverted. This condition is prevented by adding 
a constant CK to each term of the main diagonal at line 250. This 
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alters the matrix equation which is then compensated by adding 
the proper correction term to the other side of the equation. 

Assume a matrix equation of the form 

A X = B (4-32) 

It follows that 

[A + eI]X = B + ex (4-33) 

where I = identity matrix 

In the program 

DK — X 
CK = e 

Bl = B + eX line 266 

Since DK is the unknown in the matrix equation, the first trial 
value of DK for eqn. (4-33) is the previous value from the prior 
iteration. DK is then computed in Subroutine MXV which gives a 
better value of DK for eqn. (4-33) . This is repeated four times 
in lines 262 through 274, each time improving the value of DK. 

Under certain conditions CK is revised downward, lines 276 through 
280, so that 

.01 £ CK £ 1. (4-34) 

This reduces the product term eX in eqn. (4-33) so that it becomes 
less and less significant in the value of Bl. 


IV- 30 



[ 2 * 4 ] 

Subroutine MXV computes X in eqn. (4-32) without evaluating 

A For a 10 x 10 matrix this requires only about one half the 
computer CPU cycles compared with evaluating X using the 
equation 

I = A -1 B (4-35) 


MXV is modified from a generalized matrix inversion subroutine 
from McCracken. ^ 

The magnitude of adjustment vector DK can be optimized to give 

maximum reduction in the cost functional EMAX. This is accomplished 

I" 26 ] 

in Subroutine STEP using a Golden Section Search . In the 
program, the step size QS becomes a scalar multiplier on DK in 
lines 928, 952, 966, and 994. 

The boundaries of the search interval are SI and S2, with initial 
values of 0 and 1 respectively in lines 918 and 920. In the 
first part of STEP, thru line 944, S2 is increased by increments 
of 0.5, line 936, until a line thru two consecutive values of 
EMAX has a positive slope, line 934. This establishes the minimum 
value of EMAX as lying somewhere between SI and S2, as shown in 
Fig. 4-5. 

The Golden Section search is based on setting up intervals with 
ratios such that 
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(4-35) 


x + x - 1 


X = 


0 

/5 - 1 


x* - .382 


.618 


The search for minimum EMAX starts at line 948 with M = 0, 
Compute the following; 


El 

lines 948-958 
E2 

lines 962-972 


EMAX 


QS. 


EMAX 


QS 2 


. 382 (S2 - SI) 


. 618 (S2 - SI) 


(4-36) 


(4-37) 


Then compare to see which of El and E2 is larger, line 978, 


If E2 < El 


then QS 1 < QS < S2 


(4-38) 

(4-39) 


We reduce the search interval by making 


SI = QS j line 982 


(4-40) 


On the next iteration with M = 1, El = E2 and a new E2 is computed, 


E2 = EMAX j 

lines 962-972 ' QS 3 " S1 + * 618 < S2 " S1 > 


(4-41) 
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By using the ratios in the Golden Section search, it is only 
necessary to compute one new value of EMAX per iteration. After 
ten iterations, the step size QS for minimum EMAX has been 
determined. The loop exits at line 976 with M > 9. The scalar 
multiplication 

DK (I) = DK (I) *QS 

1=1,2, K (4-42) 

occurs at line 994, followed by a return to the main program. 

All the state variables should remain positive. If any elements 
become negative, they are forced positive at line 294. This 
completes the major iterative loop of Program NHMODEL. A detailed 
flow diagram for the Program NHMODEL is shown on Figures 4-6 
through 4-11, and a program list is provided in Table 4-6. 


IV-33 



NHMODEL 


QD 


i 


Li 





! ^ 


Entry 

1. Read input data NCP, NRP , NRZ , 

estimate of ZK vector, and density profile 


2. Set CK =1, QS = 1, DK vector = 0, 
and iteration count = 0 


3. Calculate K, Kl, and K2 


4 . Store the initial estimate of the ZK vector 
in PZK 



5. Compute the chanqe of variable 
S(I) * ZK (K2) - Hl'J (1,1) 



6. Compute the approximation to N(H) 
CALL CALC [AMP , HN , S , ZK] 

Return with AMP 



;L;.l 

•0 




YES / '"\ 
) 


NO 



(b) 


7. Compute the discrepancy vector, E 
CALL COST [ AMP , E , HN , ZK , EMAX , EXM] 
Return with E, EMAX, EXM 


8, Increase iteration count by 1 


9. Is EMAX < 0.01? 


10. Display iteration count, CK, EMAX, and QS 


11. Is QS < 0.01? 

12. Is CK = .01 and EMAX < .025 and QS < .03? 


Fig. 4 “6. Flow Diagram for Program NHMODEL 




NHMODEL 



13. Compute partial derivative matrix, D 
CALL GEND [D, S , ZK] 

Return with D 


T 

14. Compute B = D E 


T 

15. Compute A = D D 


16. Add CK to the main diagonal 
of the A matrix 



17. Compute approximate solution to 
DK = (D t *D + CK*I) _1 * (DT*E) 

CALL MXV t A , DK , Bl , T , K , K1 ] 

Return with DK 

18. Is EMAX > 0.1? 


19. Replace CK with 0.3 CK 


20. If CK < 0.01, set CK = 0.01 


21. Compute the multiplier, QS , to optimize 
the correction for the ZK vector. The 
new values for E and EMAX are computed. 
CALL STEP [AMP , DK , E , HN , S , ZK , EMAX , EXM , QS ] 
Return with DK, QS 

22. Replace ZK with ZK + QS*DK 



22. If any component of ZK is negative, then 
reverse the polarity of that component. 



Fig. 4-6. Flow Diagram for Program NHMODEL (Cont) 




NHMODEL 



GlD 


24. Display HN (I, 1) , HN (I , 2) ,AMP (I) ,E (I) , 
I = 1 to KL 

25. Exit 


Fig. 4-6. Flow Diagram for Program NHMODEL (Cont) 




NHMODEL 



Start 


1. Compute the effect of the complex poles 

NCP 

DEN = JJ / jPj P * \ (See Fig. 4-1) 

J-l ' J 

2. Compute the effect of the real poles 

(NCP+NRP) 

Replace DEN with DEN * JJ P-j- j 

J=NCP+1 

Compute the effects of the zeros 
AMP (L) = IS(1) ] NRZ /DEN 

4. Compute scale factor ZK(K1) = HN (10 , 2) /AMP (10) 

5. Is L = I<L? 

6. Replace AMP (L) with AMP (L) *ZK (Kl) , L = 1 to KL 

7 . Return 


8 . Increase L by 1 


Fig. 4-7. NHMODEL Flow Diagram - Subroutine 


CALC [AMC,HN,S,ZK] 






NHMODEL 



Start 

1. Compute ZKP (I) = [ZK(I) ] 2 , I * 1 to K 

2. DO 100 I - 1,K 

3. Set L = 1 

4. Compute SI = [S(I)] 2 

5. Compute S2 = S ( I ) + S(I) 

6. DO 40 J = 1,NCP 

7 . Compute Dl = [ZKP(L) + ZKP(L+1) - SI] 2 

8. Compute D2 « [ZK(L) * S2] 2 

9 . Compute RC = -Dl - D2 

10. Compute D(I # L) = 2 * ZK(L) * [ZKP(L) 

+ ZKP (L+l) + SI ] /RC 

11. Compute D (I ,L+1) * 2 * ZK(L+1) * [ZKP(L) 

+ ZKP (L+l) - S1J/RC 

12. Increment L, L = L + 2 

13. Is J * NCP in DO loop? 


Fig, 4-8. NHMODEL Flow Diagram - Subroutine 

GEND [D,S,ZK] 












NHMODEL 


60 

100 



14. Is NRP * 0 


15. DO 60 J = 1,NRP 


16. Compute RR * -SI - ZKP(L) 

17. Compute D(I,L) * ZK(L)/RR 

18. Increment L, L = L + 1 

19. Is J * NRP in DO loop 

20. Is I * KL 

21. Return with 


3 £nAMP ( 1 ) 

3£nAMP (1) 

3£nAMP (1) 

3ZK (1) 

3ZI< (2) 

3 ZK (K) 

3£nAMP (2) 

3£nAMP (2) 

.... 3 £nAMP (2) 

3ZK (1) 

3ZK (2) 

3ZK (K) 


3£nAMP (KL) 

3£nAMP (KL) 

.... 3£nAMP (KL) 

3ZK(1) 

3ZK(2) 

3ZI< (K) 


Fig. 4-8. 


NHMODEL Flow Diagram - Subroutine GEND [D, S, ZK] (Cont) 





NHMODEL 



Start 


1. Load A matrix into array T 

T(I,J) = A(I,J) 

X *" i , 2 ; • • « « / N 
J ■* X f 2 , » « « • ; N 


2. Load vector B into last column of T 
T ( I , Nl ) = B (I) , 1 = 1, 2 , , N 


3. DO 100 for K = 1 to N 

In the kth column, s tarting with the K+l 
row, determine which row of the T matrix 
has the term with the greatest magnitude. 
This is indicated as the L th row. 


4 . Is K = N ? 


5. Is L = K? 


6 . Interchange the L*" and K rows in the 
T matrix for all columns from the K th 
through the (N+l)*-* 1 , j = k to N+l 



Fig. 4-9. NHMODEL Flow Diagram - Subroutine MXV[A,X,B,T,N,Nl] 







NHMODEL 



7. Replace the T (K, J) term with T (K, J) /T (K,K) 
J * K+l to N+l 


8. Is K = 1? 


4 - 

9. For the 1^ row, from I = 1 to K-l , for 
J = K+l to N+l, replace T (I , J) with 
T (I , J) - T ( I , K) *T(K,J) 


11. For the I*" 1 row, from I = K+l to N, 

for J = K+l to N+l, replace T(I,J) with 
T(I,J) - T (I,K) *T(K, J) 



12. X(I) = T (I ,Nl) , for I ■ 1 to N 

13. Return 


Fig. 4-9 . NHMODEL Flow Diagram - Subroutine MXV [A,X,B,T,N,Nl] (Cont ) 







NHMODEL 



Start 

1. Compute E (I) = S,n [HN (1 , 2) /AMP (I) ] , I = 1 to KL 

KL 

1 ITi 

2. Compute EMEAN = E (I) 

1=1 

_ _ . „ VM (EMEAN) 

3 . Compute EXM = e 

4. Replace ZK(K1) with ZK(K1)*EXM 

5. Replace all E(I) with E(I) - EMEAN 

6. Set EMAX = Max absolute value of E(I) 

7. Is EMAX <0.5? 

8. Set ER = 0.5/EMAX 

9. Replace all E(I) with E (I) *0 . 5/EMAX 

10. Return 


Fig. 4-10. 


NHMODEL Flow Diagram - Subroutine COST 
[AMP , E , HN , ZK , EMAX, EXM] 







NHMODEL 



Start 

1, Set El ■ EMAX, M * 0, Si = 0 , S2=l 

2. Set QZK(I) = ZK (I ) , I = 1 to K 


3. Replace ZK(I) with QZK(I) + DK(I) *S2, I = 1 to K 


4. Call Subroutine CALC [AMP,HN,S,ZK] 

Return with AMP 

5. Call Subroutine COST [ AMP, E, HN, ZK, EMAX, EXM] 
Return with E, EMAX, EXM 

6. Ts EMAX > El? 


7. Increase S2 by 0.5 

8. Is S2 > 2,5? 


9. Set El = EMAX 


10. Is S2 = 1? 

11. Set SI — — 1 

12. Compute QS * SI + (S2-S1) *0.382 

13. Compute ZK (I) = QZK(I) + DK ( I ) *QS , I * 1 to 20 

Fig. 4-11. NHMODEL Flow Diagram - 


Subroutine STEP [AMP, DR, E, HN, S, ZK, EMAX, EXM, QS] 






NHMODEL 




Call Subroutine CALC [AMP,HN,S , ZK] 
Return with AMP 


Call Subroutine COST [ AMP, E,HN,ZK, EMAX, EXM] 
Return with E, EMAX, EXM 


Set El = EMAX 
Is M > 0? 


Compute QS ~ SI + (S2-S1) *0 . 6 ] 8 

Compute ZK (I) * QZK(I) + DK(I)*QS, I = 1 to K 

Call Subroutine CALC [AMP ,HN, S , ZK] 

Return with AMP 

Call Subroutine COST [AMP , E , HN , ZK , EMAX , EXM] 
Return with E, EMAX, EXM 

Set E2 = EMAX 


Increase M by 1 
Is M > 9? 

Is E2 ■> El? 

Set El = E2 

Replace SI with SI + (S2-S1) *0 . 382 
4-11. NHMODEL Flow Diagram - 


Subroutine STEP [AMP, DK, E, HN, S, ZK, EMAX, EXM, QS] (Cont) 








NHMODEL 


Page 11 


28. Set E2 * El 

29. Replace S2 with Si + (S2-S1) *0 . 618 

30. Replace DK(I) with DK(I)*QS, I = 1 to K 

31. Return 


Pig. 4-11. NHMODEL Flow Diagram - 
Subroutine STEP fAMP,DK,E, HN,S, ZK, EMAS, EXM, QS J (Cont) 
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RAM LIST 

NHMODEL 1-13-70 PAGE 1 

LIST 

hP<5 


100* 

C? 

PROGRAM NHMODEL 

10?.* 

c* 


10 4* 

c? 

PROGRAM TO INITIALIZE THE ALPHA VECTOR 

106. 

c? 

INPUT IS A SET OF N(HCD) DATA POINTS IN ARRAY HN 

10*. 

C t 


l in. 

c* 

ETA = 7KCK+2) 

no. 

C? 

HN = PRESCRIPTION VECTOR <H,N) 

114* 

C; 

K = NUMBER OF TERMS IN APPROXIMATION VECTOR ZK 

116. 

C? 

KL = NUMBER OF TERMS IN PRESCRIPTION VECTOR HN 

! ! * • 

C* 

RHO = ZKCK+l ) 

1 ?0 . 

C: 

S r ETA - HN C I , 1 ) 

1?? . 

C? 

7K = ALPHA VECTOR 

! ?4 . 

C? 


I?6. 


DIMENSION TITLE (7), ZK(1 4) , HN<20,2) , S<20) , DK( 12) 

1?*. 


DIMENSION E (20) , B(12) , A(12,12) , D(20,12), T(12,13) 

130. 


DIMENSION PZKC14), AMPC20) , W(20) , B1C12) 

13?. 

c? 


1 34, 


COMMON NCP, NRP, NRZ, K, K1 , KL 

l|6. 

c* 


13*. 


DATA NCP, NRP, NRZ /6, 0, 3/ 

1 40. 


DATA CZK(K), K = 1,14) / .06 , 7 . , .09 , 6 . 1 1 , .05 , 4 .68 , .3,4.78, . 

- 


50,3.2*,1 .02,1 .75,1 .,5./ 

1 42. 


DATA KL. HM(I,2) , 1 = 1,8) /16, 2 .09966 , .4 1 2 1 1 , 2. 

04*12, ,43*46, i *9383?, .50675, 1 . 7361 5 , .67384 , 1 .53350, .937 



1 1 , 1 .35285,1 .28677, 1 .16740,1 .90108, 1 .04249,2.64945,/ 

1 44 . 


DATA ( H N ( 1 , 1 ) , HN< 1 ,2) , 1=9, KL) / .84938 , 4 . 7 1 84 1 , .65850,9. 
*6743, .55223,1 7.92012, .48818,28.10641 , .44505,41 .50047, . 



41461 ,56.1 8123, .38580,74.68321 , .36753,92.61 606/ 

! 46 . 

c? 


14*. 


ITCNT = 0 

150. 


CK = 1 . 

152. 


QS = 1 . 

154. 


DO 70 I r 1 ,20 

1 56. 


Ed) = .01 

1 6*. 


IF (I .LE. 12) DK ( I ) = 0. 

1 60 . 

70 

CONTINUE 

162. 


K = NCP+NCP+NPP 

16 4. 


K 1 = K+l 

166. 


K2 = K+2 

i c n 
1 ® 


DISPLAY "NCP =" , NCP, "NRP =” , NRP, "NRZ =", NRZ, "ETA =", 
7K ( K2) 

170. 


DISPLAY " " 

1 72 . 


DISPLAY "ITCNT CK EMAX QS" 

17 4. 


DO I : 1 ,K2 

176. 

*0 

P7K ( I ) = ZKCI) 

1 7? . 


pr *01 = 1 , KL 

! <4° . 

on 

SCI) = 7K(K?) - HN ( I , 1 ) 

1 *° . 

C* 


1 °4 . 

i o 

C* 

o 

COMPUTE APPROXIMATION TO N ( H) 

I ■ ' ■ 0 

d*. 

CALL CALC [ AMP,HN,S,7K ] 

1*0. 

c* 


I *0 . 

c? 

COMPUTE DISCREPANCY VECTOR E 

1*4. 

c- 


1*6. 


CALL COST t AMP, E,HN,ZK, EMAX, EXM] 

1 

1 00 

CONTINUE 


Table 4-6. Program List for Program NHMODEL 
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?-■ 007 AM LIST 
t 

LIST 200:313 


200 . 


ITCNT = ITCNT + 1 



IF (EMAX ,LT. .01 ) GO TO 200 

2^4. 

C? 

DISPLAY " 7 K , ZK 

20 S. 


DISPLAY ITCNT, CK, EMAX, QS 



IF CCS .LT. . 01 ) GO TO 200 

2 70 . 


IF (CK ,EQ, .01 .AND. EMAX ,LT. .025 .AND. OS .LT. .03) 



TO ?no 

- 1 * . 

C; 


2 1 ° . 

C? 

GENERATE PARTIAL DERIVATIVE MATRIX D 

0 ! 4 . 

Cs 


01 


CALL GENU [D,S,ZKJ 

o 1 c . 

Cr 


"on * 

C* 

MULTIPLY BOTH SIDES BY D TRANSPOSE 

OO'? 

C* 


4 # 


DC 120 I = 1 , K 

99 ^ . 


PCI) = 0. 

n rl r- 


DO 120 L = 1 ,KL 

o xn m 

120 

PCI) = BCD + D(L,I)*E(L) 

?3? . 


DO 1 30 I r 1 ,K 

? 3 4 # 


D2 130 J = I ,K 

o ** <r 

v. «> V-< • 


A ( I , J) : 0. 

2 3°. 


DO 130 L = 1 ,KL 

?40 . 

1 30 

A ( I , .! ) - f, ( i , ,i > -•!•• D C L , T ) *n ( L , J) 

242. 

O 


4 4 * 

C: 

ADO CK TO MAIN DIAGONAL AND COMPLETE THE A MATRIX 

246. 

C? 


24®. 


DO 1 40 I r 1 ,K 

-? s n # 


AC J,I) r A (1,1) + CK 

9 r v,? # 


DO 140 J = 1 ,K 

£6 4 . 

1 40 

A ( J , I ) = A C I , J) 

0*36. 

C: 


25®. 

Cf 

COMPUTE ADJUSTMENT VECTOR DK ON ALPHA VECTOR ZK 

260 . 

C f 


262. 


DO 160 J = 1,4 

£64. 


DO 1 50 I = 1 ,K 

2^6. 

1 50 

Pl(I) = B(I) + DK ( I ) *CK 

2^° . 


CALL MXV I A , DK , B1 ,T,K,K1 ] 

S 70 . 

C* 

DISPLAY " M 

272. 

c* 

DISPLAY "DK =" , DK 

*7*. 

1 60 

CONTINUE 

27£. 


TR (EMAX .GT. .1 ) GO TO 1 70 

£>•>•». 


CY r CK*0 .3 

2*^0 . 


IF (OK .LT. .01) CK = .01 

0,7 4 . 

r? 


~~6 , 

c? 

COMPUTE STEP SIZE FOR ADJUSTMENT VECTOR DK AND ADD D 



K TO ALPHA VECTOR ZK 

ft 

r* 


"'r O 

« 

170 

CALL STEP r AMP, DK ,E, HN ,S, ZK , EMAX, E'XM, QS J 

, > r- «t% 


DO l on T r 1 ,K 

^ 4. 


TF ( 7 K ( T ) .IT. 0.) ZK ( I ) = -ZK(I) 

: r c . 

l on 

CONTINUE 

r" -r 

• 


no to ion 

>20 . 

200 

CONTINUE 

.30 / . 


DT SPLAY " " 

VO 4 . 


DISPLAY " J PZK ZK” 

7 r> c 
.r .« O « 


WRITE (1,021) (J, PZK(J), ZKCJ), J = 1,K2) 

$n & # 


DISPLAY " " 

2 12. 


DISPLAY " I H N ( 1 , 1 ) H N ( 1 , 2 ) AMP(I) E ( I ) " 

312. 

Table 4 
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v-/ 

UJ 




H ^ 


/--s rs 

t 
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U3 


® in 

9 9 # 



* 

«fc 

o • 

O C in 



/-S 

h* i 

— o 

__| M _* 



J—4 
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U- U- u. 



W 

|TS 

03 U <? K} CO 



CL 

o O 

s-\ 




h: 

•— « 

|p, S* ®> 





o, c 

h m it. 

in m in 

UJ 



«v 

Cl h h 

H H H 

CO 


il 

— • w ^ 

<^/ v-/ V-X 
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o 



w 



k: 


rv 

52: * f— * f— • $— 

X 


i — h 

Lti h <r 

«ac <t <- 

C <. C 

z 


w 

i— « t~* CL FI 

F.: F f: 

5l F ^ 



CL 

i-h s o tr 

r: cr C" 

Or 0: or CO 



FT 

Cl C H d 

o c o 

C C C IS: 



C 

3 o o: tr U. u.. t 

lx. La. lx. U-! 




C C 

— Co *o 

O C — 
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o; o 

c o c 

— 00 c , 

CO 

o 


CSJ o 

C V o 

o c? t 

>— j 

sr 





•4 

• * 
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K? 





O' 






cr 

{— 
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Cj 

tfl 

S3 

KO Cf c 0 

sr t 

r c S7“ v. 
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LIST 500 • 5°? 


500. 


SUBROUTINE CALC [ AMP ,HN , S , ZK ] 

50? . 

C • 


60 4. 

O 

CALCULATE NCS) OF THE APPROXIMATING FUNCTION ZK 

506 . 

Cr 


50° . 


DIMENSION AMP(20) , HN(20,2), SC20) , ZKC14) 

510. 


COMMON NCP, NRP, NRZ, K» K1 , KL 

5!?. 

O 


51 4. 

C t 

DISPLAY ” I", " 01", " DEN" 

5! ?. 


DO 50 L = 1 ,KL 

51 


1 = 1 

5?0 . 


DEN = 1 .0 

5??. 

c- 


5'-/i . 

c* 

COMPUTE AMPLITUDE RESPONSE DUE TO COMPLEX POLES 

5 If . 

cor* 

' ’ * 

c* 

DO 20 J = 1 t NCP 

5 ’ 0 . 


A 1 r ZK ( I ) 

C x o 

o * 


31 = ZK ( 1+1 ) 

534 . 


R 1 = CBI-S<L))*<B1-S(L)> 

536. 


R2 = (Rl+S(L))*(B1+S(L)) 

53". 


= Al*Ai 

540. 


Ql = SORT! ( R1+R3) *(R2+R3) ) 

5 4?. 


DEN r DEN * Q! 

544. 

C: 

DISPLAY I » Ql 9 DEN 

546 . 

?.0 

1 = 1+2 

54". 

C; 


550. 

C? 

COMPUTE AMPLITUDE RESPONSE DUE TO REAL POLES 

55? . 

C; 


55 4 . 


IF (NRP . EQ . 0) GO TO 40 

556 . 


DO 30 J = 1 , NRP 

55 c . 


A 1 = 7K(I) 

560 . 


R 1 = 5(L)*S(L) 

562. 


P2 = A1*A1 

564. 


01 = SORT t R 1 +R2 ] 

56 6. 


DEN = DEN * 01 

5 6° . 

C • 

DISPLAY I , 01 , DEN 

OO . 
> « 

30 

T r T + J 

COO 

’ ■ a 

4 0 

CONTINUE 

<•> ? /! # 


AMP(L) = S ( L) ** NR7/DEN 

57 r. 

C; 

DISPLAY I , Ql , DEN, AMP (L) , L 

C "» r* 

' • « 

60 

CONTINUE 

5 C R . 


*K(KI) = HN( 1 0 , 2 ) / A MP ( 1 0) 

5 aa 


DO 60 L = 1 , KL 

5*4. 

60 

AMP (L) = A MP(L) *ZK (K 1 ) 

S r>r 

; *. « 


RETURN 

c .*? 17 

. a 


END 


Table 4-6. Program List for Program NHMODEL (Cont) 
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LIST 600j€*<= 
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600 . 


SUBROUTINE GEND ID»S,ZK] 

60?,. 

C? 


604. 

C; 

GENERATES D MATRIX OF PARTIAL DERIVATIVES 

60 6. 

C* 


60P. 


DIMENSION SC20), 7KC14), ZKPC12), D(20,12) 

r m. 


COMMON NCP, NRP, NRZ, K, K1 , KL 

*ip. 

c* 


614. 


DO 10 I : l,K 

616. 

m 

7 KP(J) = ZK(I)*7K(I) 

61°, 


DO 100 I.; I ,KL 

^n. 


L = 1 

*0? 

6’ p 


SI r S(I)*S(I) 

6 C 4 » 


S2 = S(I) + SCI) 

6*6. 

Ci 


6 2 17 , 

C? 

TERMS FOR COMPLEX POLES 

630. 

c? 


632 . 


DO 40 J = 1 , NCP 

634. 


D 1 = CZKPCL)+ZKPCL+I)-S1)**2 

636 . 


D2 = C ZK C L) *S2) **2 

63°. 


RC r -D1-D2 

6 40 . 


D(I,L) r C?«*ZKCL) *CZKPCL)+ZKPCL+1 )+Sl > )./RC 

642 . 


DCI,L+1) : C2.*ZKCL+1)*CZKPCL)+7KPCL-H)-S1))/RC 

6 4 4 . 

40 

L = L + 2 

646. 

C? 


6 4c?. 

C: 

TERMS FOR REAL POLES 

6,50 . 

C: 


652. 


IF ( NRP »E0 . 0) GO TO 100 

654 . 


DO 60 J r 1 , NRP 

656. 


RR r -Sl-ZKPCL) 

65P. 


DCI,L) = 7K CL) /RR 

660 . 

60 

L = L + 1 

66? . 

100 

CONTINUE 

664 . 


RETURN 

666. 


END 
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P 0 1 1 A M LIST 
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T 

>1, 1ST 70p? 7<=<3 





700 


SUBROUTINE MXV [A, X, B, T, N, Nil 

n OP 


Cr SIMPLIFIED MATRIX INVERSE SOLUTION PROGRAM 

70 4 


C* GIVEN A * X = B FIND X = A**-l * B 

7Q6 


C? LOAD A C N, N) IN FIRST N COLUMNS OF T(N,N1) 

700 


C: B IN LAST COLUMN 

7 ! 0 


C: AFTER GAUSS-JORDAN PROCESS X WILL APPEAR IN LAST COL 



UMN OF T. 

7 } 9 


C; 

71 4 


DIMENSION A C 1 2 , 1 2 ) , XC12) , B(12) , T(12,13) 

*»1 6 


C; 

71 9 


DO 10 I = 1 ,N 

700 


0 0 10 J = i ,N 

r' r? 

1 0 

TO , U) : A(I, J) 



DO 20 I = 1 , N 

*7 **? f 

rn 

T ( I , N1 ) = B(I) 

7 0(7 


DO 100 K = 1 ,N 

73-0 


K 1 = K + 1 

"30 


IF (K .EQ. N) GO TO 50 

7.54 


L = K 

7 3-' 


DO 30 I r K 1 , N 

73 c 


IF (ABSCTd»K> I *GT » ABSI T (L,K) 3) L s I 

740 

30 

CONTINUE 

7.40 


IF (L .EQ. K) GO TO 50 

744 


DO 40 J r K , N 1 

7.4 4 


Q : T CK , J) 

'7 / C 


T(K,J) = T(L,J) 

"50 

40 

T(L,J) = 0 

75? 

50 

CONTINUE 

75 4 


DO 60 J = K 1 , N 1 

75,6 

60 

TCK , «!) = T(K,J)/T(K,K) 

75° 


IF CK , EQ . 1) GO TO SO 

7?Q ' 


KM1 = K - 1 

"if?. 


DO 70 I r 1 , KM1 

7 6*4 


DO 70 J r K 1 , N 1 

Iff 

70 

TCI,J) = T ( I , J) - TO ,K) * TCK, J) 

If” 


IF (K .EQ. N) GO TO 100 

770 

C-Q 

CONTINUE 



DO on 1 r K 1 , N 

- *7 /, 


DO QQ J r K 1 , N 1 

- • /' 

OO 

TCI, J) = TCI , J) - TCI ,K) * TCK, J) 

•77C 

100 

CONTINUE 

7°0 


DO HO I = 1 , N 


1 1 0 

XCI) = T C I , N 1 ) 

77/! 


RETURN 



END 
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PROGRAM LIST 
t 

LIST oQOfOocj 


NHMODEL 


1-13-70 


PAGE 7 


eg? » 


SUBROUTINE COST t A MP , E„HN , ZK , EMA X, 

-7 04 , 

Ct 




DIMENSION AMPC20) , EC20) , HN(20,2) 

eO© 

e 


COMMON NCP, NRP, NRZ, K, K1 , KL 

"10. 

Ci 




SUM = 0. 

*i 4 . 


DO 10 I : l,KL 

*i 6. 


ECI) r ALOGIHNCI ,2)/AMP(I) ] 

«i c . 

10 

SUM r SUM + ECI) 

rtf!, 


EMEAN = SUM/KL 

r-0 9 t 


EXM r EXPf EMEAN ] 

4 . 


" 7 K C K 1 ) = ZKCKl )*EXM 

rr H c 


EMAX r 0. 

©?* . 


DO 20 I s 1 , KL 

030 . 


ECI) r ECI) - EMEAN 

*32 . 


01 = ABSIECI) ] 

*34 . 


IF CEMAX ,LT. 01) EMAX = 01 

*36. 

20 

CONTINUE 

°3 * . 


IF CEMAX .LT. .5) GO TO 40 

R$0 . 


£R r .5 /EMAX 

*42 . 


DO 30 I = 1 ? KL 

*4 4 . 

30 

ECI) = E C I ) *ER 

*46,. 

40 

RETURN 

*4 r . 


END 


Table 4-6. Program List for Program NH MODEL (Cont) 
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PROGRAM LIST NHMODEL 1-13-70 PAGE 8 

t 

LIST *00?o*<? 


*02. SUBROUTINE STEP [ AMP,DK , E, HN , S ,ZK , EMAX, EXM, 8S 3 

* 04 . C« 




DIMENSION AMP<20), EC20), HNI20.2) 

-HP # 


DIMENSION D K ( 1 2 ) , QZKC12) 

*10. 


COMMON NCP, NRP, NRZ, K, K1 , KL 

* \ 2 . 

C* 


*14. 


El : EMAX 

~16. 


M = 0 

*i p-. 


SI = 0. 

-on . 


S2 r 1 . 

po? . 


DO 10 I r 1 ,K 

*24 . 

10 

P"K(I) = ZK(I) 

pgf , 

ro 

DO 30 I = 1 ,K 

-6c> # 

30 

T ( I ) = QZK(I) + DK ( I )*S2 



CALL CALC [AMP,HN,S,ZK3 

- 32 . 


CALL COST [ AMP, E,HN,ZK, EMAX, EXM3 

c34. 


IF ( (EMAX-E1 ) »GT. 0.) GO TO 40 

f ■ 3 2 . 


S2 r S2 + 0,5 

f3 n . 


IF CS2 , GT , 2,5) GO TO 100 

o40 . 


El r EM AX 

o42 . 


GO TO 20 

044. 

40 

IF (S2 ,EQ. 1 .) SI = -1 , 

*46 , 

c? 


*4*. 

100 

OS = SI + CS2~S1)*.382 

*50. 


DO 120 I = 1 ,K 

0 5? . 

1 20 

ZKCI) = QZK(I) + DK ( I ) *QS 

of 4. 


CALL CALC [ AMP,HN,S,ZK ] 

o56. 


CALL COST [AMP,E,HN,ZK,EMAX,EXM3 

o5P , 


El = EMAX 

-60. 


IF (M .GT. 0) GO TO 160 


1 30 

QS = SI + (S2-S1 )*.61 F 

06 4 , 

1 40 

DO 150 I = 1 ,K 

-66, 

1 50 

?K ( I ) = QZK(I) + DK ( I)*QS 

«-r< 7 . 


CALL CALC [AMP,HN,S,ZK 3 

p7n. 


CALL COST t AMP, E,HN,ZK, EMAX, EXM] 

* 7 ? , 


E° = EMAX 

-7 4 . 

1 60 

M r M + 1 

-76 . 


IF (M .GT. *) GO TO 200 

-70 . 


IF (E2-E1 ) 170, 180, ISO 

ocn , 

1 70 

El r E~ 

0 CO 


SI r SI + ( S2-S1 ) * .382 

-* 4 . 


GO TO 130 


1 <*0 

E8 r El 

n C ^ 

• 


S2 = SI + (S2-S1 )*.61 8 



GO TO 100 

0 r 

200 

DO 210 I = 1 ,K 

# 

8 ! 0 

C T) r DK C I ) *PS 

no r 


P8TUPN 

r* pp rp 


r M n 


S(20) , 


ZKC14) 


Tabli* 4-6. Program List for Program NHMODEL (Cont) 
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4.5 THE PSEUDOINVERSE 


f 2 8 1 

Kinkel 1 J has developed an efficient algorithm for computing 
the pseudoinverse of a matrix which can simplify the solution 
of eqn. (4-32) . 

The weighted least squares solution for X using the pseudoinverse 
is represented as follows: 

X = [A T WA] + A T WB (4-42) 

Even if A is poorly conditioned so that it cannot be inverted, 
the weighted least squares solution of X is given by eqn. (4-42) . 

The use of the pseudoinverse will avoid the iterative method 
associated with eqn. (4-33) . 
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SECTION V 


HORIZONTAL PROCESSING 


5 . 0 INTRODUCTION 

Present methods for routine topside ionogram data reduction 
depend on the following simplifying assumptions: 

a. Vertical propagation of the echo pulse 

b. Spherically stratified ionosphere over the period of the 
X-Trace 

It is a basic property of the topside sounding satellites using 
swept frequency radar techniques that successive h'(f) echo 
pulses correspond to ionospheric soundings at different earth- 

based coordinates. It is well known that the ionosphere is 

f 1 4 ] . 

nonspherically stratified and since the Alouette II satellite 

moves at an average velocity of approximately 7 km/sec the 

effects of horizontal gradients in the orbital plane should be 

taken into consideration. 

Present ionogram data reduction methods for h' (f) — *> N(h) 
computation assume a constant electron density profile for the 
ionosphere during the period of useful data acquisition in each 
ionogram. In addition, the Forward, Inverse, and Matrix processing 
methods are concerned with reducing the data from only one 
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ionogram at a time. In the Forward and Inverse methods, data 
reduction begins at satellite height and works downward as 
laminations are added to make up the N(h) profile. This has 
become known as vertical processing. 

The method of Horizontal Processing described in this section 
is proposed as a method of providing first order correction of 
the effects of horizontal gradients. The method assumes 
a constant horizontal gradient between adjacent ionograms 
in the same pass at equal true height lamination boundaries. 


In this section the horizontal gradient is defined as the rate 
of change of electron density with respect to distance in the 
orbital plane at constant true height. 


This can be represented as follows: 



3 

where N = density in electrons/cm 

x = distance in orbital plane 

3 

h = true height in 10 km 


( 5 - 1 ) 
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Since time is one of the parameters associated with an ionogram, 
it is more convenient in data reduction to work with the gradient 
as a function of time than with distance. The gradient with 
respect to time is as follows: 

(H) h - (f ) h (n) h 

At any true height level, the term (^) is essentially constant 

h 

for the period of one ionogram. For Horizontal Processing the 
gradient in eqn, (5-2) is assumed constant between adjacent 
ionograms at each lamination boundary. 

Representative N(h,t) profiles for two ionograms are shown in 
Fig. 5-1 referred to a time scale. The N(h) profile shown as 
curve AB is the intersection of the ionospheric density surface 
ABGH and constant time plane t.^. The N(h,t) profile of Ionogram 1 
is shown as curve AC which is the intersection of ABGH and a 
skewed time surface. This is the type of electron density profile 
that will be computed with the Horizontal Processing method. 
Similarly the N(h,t) profile of Ionogram 2 is shown as curve EG. 

With Horizontal Processing it is necessary to work with a minimum 
of two adjacent ionograms from the same pass. The horizontal 
gradient is then included in the N(h,t) computation by linear 
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interpolation of N with respect to time along each lamination 

boundary. This condition requires that the lamination boundaries 

be at the same time height for each ionogram. Horizontal Processing 

is most easily implemented by using Inverse Processing because the 

heights of the lamination boundaries will be the same for all 

ionograms in a pass, except at satellite height. We can also 

assume that 9N/3t is a constant along the satellite track so that 

N(h ,t) can be obtained by linear interpolation, 
s 

5.1 HORIZONTAL INVERSE PROCESSING 

The method of Horizontal Processing will be developed as an 
extension of Inverse Processing, which is described in Section III. 
X-Trace data from two adjacent ionograms is scaled in conventional 
manner, and then curve fit so that h'(f) is defined for any 
frequency on both X-Traces . 

The value of N(h^,t) at satellite height is computed for each 
ionogram using the relation 

f N - / f )«i - 'll' !V :i 

where fj. = f xs 

then N(h 1 ,t) - 12400 f^ 2 (5-4) 
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Referring to Fig. 5-1* the electron density at the satellite 


for Ionogram 1 is denoted by N(h^,t^) and for Ionogram 2 by 
N (h 1 , t 2 ^) . The value of N (h^, t) at any time between t^ and 
is then determined by linear interpolation. The gradient along 
the satellite track is approximated by 


c _ AN 
1 At 


N(hi,t 2 i) - N(h 1 ,t 11 


(5-5) 


'21 


'll 


Then at any time t along the satellite track from t^ 


to t. 
16 


N(h x ,t) 


N (h. 


' t ll ) 


s 1 (t - 


til) 


(5-6) 


It is convenient to denote satellite height by h^ and yet, 
because of the elliptical satellite orbit, h^ will change from 
one ionogram to the next. This doesn’t change the validity of 
eqn. (5-6), but the change in height should be taken into 
account when computing the gyrof requency f^ for successive 
laminations . 


The bottom of lamination 1 is denoted by h ^ in Fig* 5-1. 
Ionogram 1 the first trial value of fcKh^/t-^) 9^ ven 
(3-16) which is repeated for convenience. 

b?) 

N 2 = N l e 


For 
eqn . 


(3-16) 
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where 


N x - N(h ] _,t 11 ) 


a 2 is derived empirically 


From this a value h^ ( f 21 ) as 3-2 is computed by Reverse 

Processing. We now have a trial value for t ^ 2 i n Fig. 5-1 given 
by 


'12 


t (f 21 ) 


(5-7) 


A value for N(h^,t 12 ) now com P ute< 3 with eqn. (5-6), using 
t = t-^ 2 * For t ^ le seconc ^ iteration choose N(h 2 ,t^ 2 ) using eqn, 
(3-18) 

Hj2 = Njl (l + a) 


(3-18) 


and compute a second point h^(f 22 ) as -*- n Fig. 3-2. A straight 
line thru points h'(f 01 ) and h'(f 00 ) will intersect the X-Trace 
at frequency £ 23 * Next, update the estimate for t ^ 2 by 


'12 


t( f 23 ) 


(5-8) 


The value of N(h^,t^ 2 ) is again updated using t = t 
(5-6) . For the third iteration, 


in eqn. 


N2 


/ f 23 (f 23 


f H2 J 


(3-19) 


and a third point h^(f 22 ) is computed. A segment of a circle 
passing thru the three computed values of h^(f 2 ) will intersect 
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Then as described on page III-9, 


the X-Trace at frequency f^. 
the final value of N(h 2 ,t^ 2 ) can be computed from eqn. (3-19), 

(3-20) with fj = f 2 4 * 

The corresponding point N(h 2 ,t 22 ) is confuted for Ionogram 2 
in a similar manner except that NOi^t^) is computed by 
extrapolation. Assume the same slope as given by eqn. (5-5) , 
then 

N(h l /t 22 ) = N(h l' t 21 ) + S l (t 22 " t 21 ) (5-9) 

The horizontal gradient for all other laminations is approximated 
by 


s - "‘W - H(1 W 

j ‘ 


(5-10 


'2j 




where j = 2 , 3 


, J lamination boundary number 


In computing the values of N(h.,tj. .) and N(h.,t~.) for succeeding 

3 j J ^ J 

laminations, the significant difference in Horizontal Processing 
is that the delay in previous laminations changes with the 
interpolated values of N(h^,t^) and the extrapolated values of 
N(h^,t 2 ^) where i = 1 , 2, . . . . , j. In Vertical Processing it is 
assumed that the N(h) profile from previous laminations remains 
constant. From the example shown in Fig. 5-1, it can be seen 



that correction for the effects of horizontal gradients requires 
the use of extrapolated data in Ionogram 2, while for Ionogram 1 
the correction is obtained by interpolation. The use of 
extrapolated data values can be minimized by processing the 
ionograms in a pass starting with the last two ionograms and 
working in reverse sequence. After the last two ionograms have 
been reduced this way, the one referred to as Ionogram 1 becomes 
the interpolation reference for the next previous ionogram and no 
further extrapolation is required for the remaining ionograms in 
the pass . 


The final output N(h,t) profile is then properly the intersection 
of the ionospheric density surface ABGH of Fig. 5-1 and a constant 
time plane. It then is a matter of choosing the time plane in 
each ionogram at which to compute N(h,t) . In Fig. 5-1, the curve 


AB is the ISt(h,t) profile where t^ 
select the time plane at t = t^g. 


= t ( f ) . One may wish to 
xs 

This is convenient in that the 


curve DC is already known at the time the value of is computed 
at the lower boundary of the last lamination. 
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5.2 HORIZONTAL MATRIX PROCESSING 


Horizontal processing using the Matrix method of computing the 
N(h,t) profile can be implemented by setting up the analytic 
model of N (h, t) to converge to curve AC of Fig. 5-1. This will 
then provide a direct mapping relationship as shown in Fig. 4-2, 
since the variables in both planes are functions of time. 

As a starting point, it may be desirable to compute curves AC 
and EG of Fig. 5-1 using the method of Section 5.1. Then make 
adjustments on AC using the matrix method of Section IV with 
curve EG as an interpolation reference for accommodating the 
horizontal gradient. When a matrix solution of curve AC is 
obtained, it then becomes the interpolation reference for a matrix 
solution of curve EG. 

A different starting point would be to compute N(h) for ionograms 
1 and 2 using the Inverse Mixed-Mode Processing method of Section 3.3. 
The computation for Horizontal Matrix Processing would then continue 
as described above. 

Interpolation of N(h,t) for horizontal gradient compensation 
would proceed as in Section 5.1, where reference values for 
interpolation along curve AC would be computed with respect to the 
analytic model along lamination boundaries at true height levels 
computed as in Section 4.3 
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Abstract 


This paper proposes a method which may be used to find an 
electron density profile yielding the least squares error fit 
to the observations. Operato/: judgment of the quality of 
observation is taken into account in the data reduction 
process. A method by which time and position variations 
occurring over the observation period may be corrected for 
is also proposed. 
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Introduction 


Some basic techniques for reduction of ionogram data have 

• [i 3 1 

been described by Jackson 1 J „ Improved processing tech- 
niques have been investigated by Madsen 1 J . The basic 
problem is to find the electron density, N, as a function 
of altitude, h, corresponding to the observation of apparent 
path length, h’, as a function of frequency, f. The relation 
ship between N(h) and h ' ( f ) is such that it is convenient to 
map a set of points from the N(h) plane to the h ’ ( f ) plane 
and inconvenient to map from h' (f) to N (h) . For this reason 
iterative procedures have been used to find N (h) given h 1 ( f ) • 
Madsen terms the mapping from N (h) to h ! (f) the Inverse 
mapping. Madsen * s inverse mapping technique requires inter- 
polation between observation points in the h'(f) plane . 
Neither Jackson ' s nor Madsen ’ s techniques consider the time/ 
position variation of the observation over the duration of a 
single ionogram. 

Madsen has suggested** that present methods of processing 
ionogram data might, be improved by a method incorporating 
either or both of the following features : 


* Numbers in parenthesis refer to bibliography , 

** Verbal discussion. 
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1. Time/position variation over the observation 
period of one ionogram is accounted for. 

2. Data from usable portions of all three 

traces, Z, 0, and X, from a single ionogram are used 
effectively to generate the N(h) profile. 

The balance of this paper describes proposed methods for 
accomplishing these goals. The methods described for 
accomplishing items 1 and 2 above are essentially independent; 
i.e., either may bo used without the other. 
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Approximations of h*(f,t) 

A basic problem with swept frequency ionogram data is that 
h' (f) is also a function of time (and hence position) whereas 
what we would like is h* (f) data for a fixed time and position . 
Therefore, we should write our observation as h* (f,t) rather 
than h' (f) to show the time dependence. This is illustrated 
in Fig, TR-1-1 in which we see that h' (f , t) forms a surface in 
the three dimensional cartesian coordinate system (h ' , f , t) . 

All observations made using a swept frequency sounding system 
must lie in the observation surface, f(t). Therefore, as shown 
in Fig. TR-1-1, a swept frequency ionogram consists of the locus 
of points forming the intersection of the two surfaces h 1 (f , t) 
and f(t). A sequence of ionograms can be visualized by a 
sequence of observation surfaces spaced along the t axis in 
Fig. TR-1-1. 

While it is true that the time/position dependence of the 
electron density profile, N(h,t) , could possibly be found 
directly from h'(f,t), it appears far simpler to first find 
h 1 (f,t, ) , where t, is some particular point in time, and from 
this find N (h, t^) . With this approach the mapping (or inverse 
mapping) problem is not increased in complexity by the added 
dimension of time/position. 



h 



:cy 


OBSEPVATIOl 
SUP.FACE 
f = f(t) 
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Basically/, our problem, then, is to construct the surface 
h' (f,t) of Fig, TR-1-1 based on our h ' (f (t) ) observation. The 

intersection of this h 8 (f , t) surface with the plane, t = t^ , 
yields h ! (f , t^) , the desired result. 

A curve may be fit to a set of points by a number of common 
methods including Lagrange polynomial interpolation v ' , 
spline fit ^ 3 ^ , and an approximating function ^ 3 ^ . Both the 

Lagrange polynomial interpolation and the spline fit pass 
exactly through all of the given points whereas the approxi- 
mating function will not. Therefore, for a large number of 
noisy observations the approximating function method may give 

( 4 ) 

superior results. Birkhoff and Garabedian have described 

a method of spline fit for a surface. However, because of 
the large number of noisy observations within an ionogram, 

I believe that an approximating function method would be best. 
Actually, a three stage hybrid method using an approximation 
function in two stages and a spline fit in the other will be 
described. 

In a typical one dimensional approximation problem we select 
a function g(a^,f) with which we wish to approximate the 
function h 1 (f) . The parameters, a^ are selected to minimize 
the error as defined by some measure of error E{g(a^,f), h ’ (f ) } . 
In this problem we shall account for the time dependence by 
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using time dependent parameters, a^(t). The time dependence 
of a^(t) can be found using points from sets of ionograms . 

The choice of an approximating function for g(a^ f f) is 
critical as regards the number of parameters, a^, required 
for a suitably good approximation. A set of orthonormal 
functions is frequently selected as a set of basis vectors 
for analytical convenience in evaluating a^ . However, in 
our case we shall use a gradient method ^ for finding a^ 
so that we may select a g(a^,f) so as to minimize the number 
of a^ instead. To promote convergence of our gradient method 
of selecting a^ we shall use a least squares measure of error 
(&2 norm) which should be satisfactory from a philosophical 
standpoint. 

In selecting the approximating function g(a^,f) we consider 
the general shape of an ionogram shape as shown in Fig. TR-1-2. 
We notice the following features? 


1. h # (f 1 ) = 0 


2 . 


dh' 

df 


f=f 


1 


• — oo 


3. h " ( f 2 ) 


oo 
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A typical orthonormal function expansion would require many 
terms to approximate these features. Therefore we consider 
the following functions. 

a 2 

gfa^f) = a 1 (&n f/f 1 ) , 0 < a 2 < 1, ^ < f < a 3 

g (a^ , f ) = a^ + a 5 (in f) + (In f ) 2 , a g < f < a^ 

~ a g 

g (a^ , f ) = a g (£n f 2 / f ) ? a 7 < f < f 2 


Notice that the points of the piecewise function are included 
in the parameter set. This function clearly satisfies all 
three of the features of h s (f) listed above. If we constrain 
g(a^,f) and its first derivative to be continuous at the 
joints then we will have only five parameters to adjust for 
the least mean square error fit; i.e, f 



g(a i ,a^) = g(a i ,a+) 
lf (a i^ a 7 ) = |§ <a i' a 7 ) 
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and we find { a^} such as to minimize 

[h 5 (f) - g(a i? f)] 2 

subject to the constraints on continuity* Using the four 
continuity constraints we eliminate four of the a^ 's leaving 
five a^'s which we can renumber and assign as elements of a 
minimizing vector to be found using a gradient technique. 


The steps involved in finding an approximation ares 


Find a function, g(a^j,f), approximating the j 
ionogram h ' ^ (f , t) . 


2 . Select a set of points , {g (a. . , f ) , f } lying in the 

1] IB Itl 

rjL 

j observation and apply spline interpolation in time 
between successive ionograms to find a set of curves 
{g’(f , t) } . The intersection of these curves with a 
plane of constant time, t = t^, gives a set of points, 
{ g ' ( f , t^) ) . This set of points approximates ionogram 
data obtained at the same instant in time . 


Find a function g (eu (t k ) , f ) giving the least mean 
square error fit to the point set { g ' ( f m , t^) } . This 
function g (a. (t, ) , f) is our approximation to an 

X K. 


3 . 
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ionogram taken all at the same time instant, t^. 

From this we can select a set of points, { g (a^ (t^) , £ m ) } , 
from which to compute N(h,t^) . 


In finding the approximating function in step 1 above we may 
use a weighted least squares cost function to discriminate 
between good and bad data points from the original ionogram 

+"U 

data. We define the cost function for the j ionogram data 
as 


J . 


3 



g(a. ,,f ) 

^ ij n 



where the weighting function, w^ , is made larger for well 
defined points, h'(f ), and smaller for poorly defined points. 
Assignment of w n would be by operator judgement using some 
arbitrary standard scale. 
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Approximation of N(h,t) 

The basic problem of determining N (h) from h ' (f ) data (or 
from g( ai (t k ) , f ) which approximates h ' ( f , t^ ) ) appears to be 
a problem in nonlinear estimation since we are given a set 
of noisy observations and some nonlinear mapping function 
relations N (h) and h ' ( f ) • If we were able to define a 
suitable approximating function for N (h) , we should be able 
to find a least squares error fit to the observations using 
a gradient technique. 

Examination of Fig. 3 shown in ref erence ^ 1 \ shows that a suitable 
approximating function might be 

-b-,£n N 

h = b, + b 2 e 

~b 3 

= b 1 + b 2 N 
or 



where b^ , b 2 , and b^ are parameters to be determined. These 
parameters can be written as elements of a state vector 


Z (b]_ '^2 '^ 3 ) 
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The above approximating function for N (h) can be used to 
find (using inverse mapping) h’ (f) , the h'(f) data computed 
from N (h) . A cost functional is now defined as 


N 


J “ X V h V f n> - h V f n>> 


i=l 


M 


+ > B. (h ’ (f ) - h' (f ) ) 

1 o n' o n 

i=l 


L 



i=l 


h' (f ) ) 
x n 


2 


where the subscripted z , o, and x refer to the three major 
types of waves producing the ionogram and the weighting factors 
, and would account for the relative quality of the 
data points . 


To apply the time position correction described in the 

previous section wo would replace h'(f ) by g(a^(t^),f ) in 

the expression for J. In this case the point weighting has 

already been accounted for in determining g(a. ( t, ) , f ) so 

iK n 

that A^ = in the expression for J. 
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A gradient technique v ' can now be used to determine the 
state vector 

Z = (b 1 ,b 2 ,b 3 ) 

such as to minimize the cost function, J. Therefore, we 
will have determined a least mean square error fit of our 
approximating function, N (h) , to the observations, h'(f) 
or N(h,tjJ using g (a^ ( t^) , f ) . 


! 
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Conclusion 


Methods have been described for finding N(h) or N (h,t^) from 
ionograms using approximation theory. A gradient technique 
is proposed to find the parameters of the approximating 
function. The success of the technique will depend upon the 
following factors: 

1. The appropriate selection of the approximating model. 

2. The speed of convergence of the gradient method 
computation of approximation model parameters. 

The approximation model proposed for N (h) is quite simple and 
when used with h ' (f ) data would lead to results which could 
be compared directly with present N(h) results. Therefore, I 
recommend that this portion of the proposed method be evaluated 
first. Since the proposed model is so simple an additional 
term involving two additional parameters would not impose a 
computational hardship if experiments indicate a more accurate 
model is required. 

Approximating h ’ ( f , t^ ) by g(a^(t^),f) appears to be a more 
difficult problem since estimation in two dimensions, frequency 
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and time, is required. Also 'the form of the approximating 
function appears to be more complex than for N (h) . Another 
factor is that results obtained by accounting for the 
time/displacement changes over a single ionogram could not 
be compared directly with results obtained by current 
methods. For these reasons I recommend that the approxi - 


mation of h ' ( f , t^ ) 


be attempted as a second step following 


the approximation of N(h). 
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A WEIGHTED LEAST SQUARES 
APPROXIMATION METHOD* 


Given a function f (x) in terms of a set of data points {x^, f(X^)}, 

i = 1,....,N, we wish to approximate f(x) by the function 

T 

F(X, a), a = (a^,....,a ) , where a is a parameter set which we 

will select to minimize the weighted mean square error of the 

approximation. We approximate the differential of F(X, a) at 

X. as 
1 


IF (X . , a) 


AF(X i , a) 


1 — 


9a. 


Aa^ + .... + 


3F(X. , a) 


i — 


3 a 


A a. 


K 


K 


~ V F (X . , a) Aa 
a ' l ' — ' — 


l — 1 p . « « . / N 


where V a F(X^, a) = 


3F(X ±f a) 


3F(X i# a) 


9a, 


9 • • • M 


3 a 


K 


A a = (Aa 


4 v t 


We \^ould like to find the AF such that 

f(X.j) - F(X i , a) - AFU^ a) = 0, i = 1, ...., N 

or AF(X., a) = f (X i ) - F(X. , a), i = 1, N 


Now we wish to find the Aa which will produce this result. We 
have 


*A. J. Mai linckrodt , "A Weighted Least-Squares Adjustment for Network 

Synthesis" Astrodata Technical Report, July 1964. 
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j f (X 1 ) - FCX^ a) 

f (X N ) - F(X n , a) 


e = A Aa 



i 




f(X x ) - F(X X , a) j 


Va F(X 1? a) 

. 

e = 

6 j 

* 

j © 

, A = 

| 

© 


f(X N ) - F(X n , a) 


Va F(X n , a) 


i 


L_ 


and the weighted least mean square estimate for Aa is 

Act* = (A T QA)“ 1 A T Qe 


, Va F(X 1 , a) 


S Aa 


j^Va F(X n , a) j 


where Q is a symmetric positive definite N x N matrix . Now we 
iterate on a producing sequentially 


m+1 

a 


a m + Aa* 


We continue iterating on a until 


A a* 


0 


Now we inquire as to the error criteria imposed on e by this 
method. Wo note that: 


Aa* 


0 < > A T Qe = 0 
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3 F(X , a) j* 

, X) + . 

j=l 


3 F ( X jj » 2L) 


N 


3 a, 


3 a. 


k 


X Q„ -e . = 0 

N j j 

j = l 

k - 1 ; • • • • / K 


Now consider the cost functional defined by 


J = 


N N 

£V= E E 

i=l j=l 


e . 
ij 3 


N 


N 


- e . 


E 


+ e. T E 0 xt . e . 

N ' N j j 


j = l 


j = l 


min J 
a 


(V J) 
a 


T 


0 < > 


3 J 
3 a, 


0 , k = 1 , . ... , I< 


3 J 


3 a 


k 


3e . 

1 


N N /„ 

„ / 3e . 

E E U - 1 Q- -e. + e.Q. . y- 
\^ a k - 1 i 13 3 a 

i=l j=l x 


N N 


3e . 


2 E X t? — — Q - -e. since Q. . = Q . . 

*-• ^ 3 a, 13 3 13 31 

i=l j=l 


But 


3e . 
x 


3 F ( X . , a) 

3 a, 


N N 


. 3 J 


3 a 


k 


3 F ( X . , a) 

•2 E X - 1 “ 


i=l j=l 


3a k i j j 


Q ■ > e • / k If . ... f K 
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We see that 


(VaJ) =0 


N N . 

8 F (X . , a) 

X 

Z^j JL*j 

i=l j=l 


9a k j e j 0 ' k lf * ***' K 


and since 3a* = 0 < 


N N ~„, v . 

3F(X.j_, a) 

53 Z) 3 ^ j e j 

i=l j=l 


= 0 


we have Aa* = 0 


(V J) A = 0 --=> min J 

2 . a 


Therefore, by iterating on a until we have achieved A a* = 0 we 
have actually achieved a weighted least squares fit to f(x) by 
F (X, a). The desired result is achieved by sneaking in through 
the back door . 

We see also that if J were quadratic in ot then F would be linear 
in a so that AF = VaF ° Aa would be exactly true and convergence 
to Aa = 0 would be achieved in exactly one step as with the 
Newton-Raphson method . For the case in which F is linear in a, 

J quadratic in a, the Mallinckrodt* and Newton-Raphson methods 
are identical . An advantage of the Mallinckrodt algorithm over 
the Newton-Raphson is that only first order partial derivative 
are required while second order partials are required for the 
Newton-Raphson. See Fig. TR-2-1 for a flow diagram for the least 
squares determination of N(h). 









o* J 



START : 


1. Initialize a., -j 

3 


1 ,K 
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2. Compute h. , i = 1,N, ) f, = f, T 

3. Compute i = 1,N (Madsen Inverse Mapping) 


Do 20 j = 1,K 

Perturb a . 

3 

Do 10 i = 1 , N 

Compute y 

Compute hj_j* 

10 Compute 

20 Restore a , 
3 


= f 


t 


r-* 

h . , 

13 


(Madsen Inverse 
h 


Compute 

A. . 

il 


5. Compute Aa = (A T QA) 1 A T Q(h'-h') 

6 . | Aa | < g? 

7. a = a + Aa 


8. N = F (h, a) 


TR 2 1. Flow Diagram - Least Squares Determination of N(h) 


Technical Report 3 
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INVERSE MAPPING TO SPECIFIED f 


We would like to perform the inverse mapping operation such that 
the frequency, f , of the end point is pre specified, as shown in 
Fig. TR-3-1 . 


h' (f) N (h) 



Fig. TR-3-1 


To do this we have the following relations for the X-Trace: 

f = ~ f„ + ~ / f, , 2 + 4f 2 
x 2 II 2 / H N 

Subtracting 1/2 f from both sides and squaring gives 


f 


2 


x 


f H f x 



= 0 


( 1 ) 


(2) 


According to Jackson 


[ 13 ] 


f H = 2.8 B 


N 


/ 


N 


12400 


(3) 

(4) 
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From our approximation formula for N we have 

N = F(h,a) (5) 

A representative curve of this function is shown in Fig. TR-3-2. 

N 

A 

h 

Fig . TR- 3-2 



We see that for any specified value of h we may compute 

N from eqn. (5) 
f^ from eqn . (4 ) 

f 1 7 1 

B from FIELDG 1 1 and interpolation 

from eqn. (3) 
f from eqn. (1) 

Since N and B are monotone decreasing in h , f^ and f H are also 
monotone decreasing. Therefore we conclude that f is likewise 
monotone decreasing in h . This property permits us to use a 
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successive approximation technique for finding h 3 



0 . 


A convenient algorithm for this purpose is Bolzano's root finding 
technique . 


We may elect to select h to within some given tolerance or to 


iterate until 0 < 


< e . Figures TR-3- 3 , TR-3-4 , and 


TR-3-5 show h selected to wi thin . 01% of a maximum value . 


* 

Wilde, Douglas J. , Op t imum Seeking Methods, Prentice Hall, 1964. 
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START; 



1. Initialize parameters 
*> = A = h/2 

2. Compute N , f , B, f^, f x , A = A/2 



5 . h = h ~ A 


6 . h = h + A 


7. N = F (h f ot) 
f = /N712M0 

B = from existing subroutine 
f = 2.8 B 

f = I f + — /f^ + 4f^ 

x 2 H 2 / r H N 


Fig. TR-3- 3 . Flow Diagram - Selection of h 3 (N,h) — ♦- (h; f ) 

X 

Where f is Specified 
x r 
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SIMPLIFICATION OF MATRIX INVERSION PROBLEM 


We are given the problem 


AX = b 

from which we wish to find 

X = A -1 b 

A Gauss-Jordan inversion process starts with 

(a : i : b) 

• • 

and through a series of row operations becomes 

(i : a -1 : x) 

• * 

Ordinarily the last column of the augmented matrix is omitted 
and X is found by a subsequent multiplication. 

X = A -1 b 

For the special case in which A is a symmetric matrix we may 
simplify the procedure to reduce only the upper right triangle 
of the left matrix to 0's which will produce only the upper 
right triangle and diagonal of A . We know that if A = A , 
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“1 T -1 -1 

then (A ) = A so that the lower left triangle of A is 

found from 


(A 1 ) i j = (A -1 ) ji 


• 2 
This would reduce the number of operations from N (N + 1) 

2 

for a savings of 2N - 3N. 



1 

2 


N(N + 1) 


However, when A ^ is not required for some other purpose than 
finding X, an even greater economy is available by deleting the 
I matrix augmentation, A further advantage is that A need not 
be symmetric. Here we would have 


(A : b) 


(I 


X) 


Not only would we save + 1) = ■j N (N 2 - 1) operations in the 

. . 2 

inversion process but also another N operations in the multi- 
plication, X = A ■'"b, which has been eliminated. Therefore, the 
total savings in this method amounts to N (N 2 - 1) + N 2 | 

as compared with the standard inversion and subsequent 
multiplication . 


For example, consider a 10 * 10 matrix; i»e., N = 10. The number 
of operations in finding the inverse is 


N 2 (N + 1) = 1100 
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. 2 

The subsequent multiplication requires an additional N = 100 
for a total of 1200 operations* 

. . -1 12 

Finding X directly and not A requires j N (N + 1) = 605 

operations for a savings of 595* 

+ N 2 } 

> 10 , 

Subroutine MXV for computing x, where 

Ax = b 

is listed on the following page. 


* $12 
This agrees with the expression for savings of w N(N - 1) 

990 + 10 0 = 59 5, Therefore, for N reasonably large, say N 

the savings amounts to almost 50% „ 


3 
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TIME SKEW PROBLEM 


Fig. TR-5-1 shows the basic assumption made in the Jackson method 
of ionogram reduction. Perhaps instead of assuming successive 
values of one could interpolate between successive ionograms. 
Since (electron density at satellite height) is the only value 
obtained without assumption of previous values, this would apply 
only for essentially constant satellite height between adjacent 
ionograms. Also, while interpolation would be possible in reducing 
the first ionogram, extrapolation would be required for the 
second. Interpolated and extrapolated values would be used in 
place of the assumed values of shown in Fig. TR-5-1. 


t. 


N, 


. - o ASSUMED „ N„ ASSUMED 

<j T U 

! h» (f,) 

j -L 

N x l 1_N 1 assumed _ 

h' (f 2 ) 

N 2 1 ' 


N. 


o_N q assumed 

ASSUMED.. 

„ N 9 ASSUMED 
h' (f 3 ) 


Fig. TR-5-1 
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For the least squares estimation approach an entirely equivalent 
set of assumptions is made as shown in Fig. TR-5-2 . 


t 


0 


f 


0 


O- 


t 


3 


~ f Q ASSUMED 


h' (f. ) 

V 1 


O-- 


h ! (f 9 ) 

i _ 


„ h , (f 1 ) ASSUMED 
h ’ ( f 2 ) ASSUMED 
1 / h* (f 3 ) 


Fig. TR-5-2 


Only four data points are shown for convenience. 


The concept of interpolating h'(f i ) data points from successive 
ionograrns is immediately obvious from the last figure . Unlike 
the case for values , only interpolation (no extrapolation) is 
required for h'(f^) values. 

Since the approach of interpolating h ' (f ) data is regarded with 
suspicion by knowledgeable people skilled in the art, and since 
solution by interpolation of N(h) is fraught with difficulty , 
the following approach is proposed. 
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1. Interpolate h ' ( f ) as described above. 

2. Compute N (h) from interpolated h ' (f ) data. 

3. Interpolate N(h) as obtained in step 2 in order to 
compute h ' (f ) . 


4. If h' (f) agrees with h 1 (f) , then we have found N(h) 
profiles which when interpolated yield results 
matching the physical measurements. 
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OPTIMAL STEP SIZE 
BY ONE DIMENSIONAL SEARCH 


The computation of step size using the Hessian matrix is based 
on the assumption that the cost functional can be reasonably 
well approximated by a truncated Taylor series having terms of 
no higher order than second. This assumption is frequently 
unjustified resulting in a step size which may be either too 
large or too small. For the case in which the step size is too 
small the cost reduction at each iteration is less than could 
be achieved resulting in the need for more iterations than 
necessary for convergence. A more serious defect occurs for the 
case in which the step size is too large; in this case the cost 
may increase for some iterations, seriously impairing, if not 
completely preventing, convergence to the optimal solution. 

This effect lias been noted in practice. 

Since we know that there exists some step size in the gradient 
direction which if sufficiently small will produce a decrease in 
cost, we should be able to produce a cost reduction in every 
iteration until the optimal solution is achieved. What we would 
like to do is to find the step size that will yield the greatest 
reduction in cost for each iteration. This optimal step size 
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may be found directly (as opposed to indirectly as before) by 
means of a one dimensional search technique. Before we can apply 
one dimensional search over an interval we need to define the 
interval, S*, to be searched, as shown in Fig. TR-6-1. 



> S 


Fig. TR-6-1 


For the one dimensional search technique to be effective, J(S), 

0 < S £ S* must be unimodal. Therefore, we would like to define 


S* as 


min 

S 



0 


S* 


Since we know there exists 6 > 0 such that J(<S) < J(0), it is 
clear, that defining S* as shown guarantees that J(S) will be 
unimodal oerr the interval 0 < S < S*. 


Having defined the interval, we can now apply either a Fibonacci 
or Golden Section search method.* 


*Wildo, Douglas J., Optimum Seeking Methods , Prentice Hall, 1964 
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We can approximate S* by the following algorithm: 

J(S i ) = J(z - S £ v z J) , S Q = 0 

j(s i ) < J(s i _ 1 ) ==> S* > s i , s 1 > 0 

If J(S^) < J ( S i- 1 > ' then let S^ + ^ = KS^, K > 1, and repeat test. 

Continue this process until J(S i ) J(S i _ 1 )i Then let S* = s.. 

-3 

Perhaps = 10 and K = 2 would be good choices initially. 
Having now determined a value for S* we may search this interval 
for the value of 0 < S < S* whic& minimizes J. 

A search by Golden Section will be almost as efficient as a 
Fibonacci search and should be less subject to round-off error 
in the final steps. Twelve iterations will define S to within 
0.5% of S*. 

We can compare the burden of computation for the Golden Section 
search relative to the Hessian matrix method as follows. The 
cost must be computed once per iteration for the search method 
and once for each element of the Hessian matrix. The number of 
elements, K, to be computed in the Hessian matrix is given by 

K = -|(N + 1) N 

where N is the number of state variables. 
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N 

K 

2 

3 

3 

6 

4 

10 

5 

15 


If we assume that twelve iterations of the Golden Section search 
is adequate, then we see that for four or less state variables, 
the evaluation of the Hessian requires fewer computations but 
more computation for five or more state variables. 


Refer to Fig. TR-6-2 and Fig. TR-6- 3 for flow diagrams for the 
interval computation and the Golden Section Search for optimal 
step size. 


Fig . TR-6- 4 is used in conjunction with Fig, TR-6-1 to show the 
sequence of computing values in successive iterations in the 
Golden Section Search . 
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START: 



E2 « El 


2. Z i = Y i " SD i' 1 = 1/ , N 

Compute E 

3. E > E2? 

4. S = 2S 
E2 = E 

5. S B J(S i ) > J(S i _ 1 ) --■> > 0 


Fig, TR-6- 2 . 


Flow Diagram - Interval Computation 
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START : 

Y 1# S 

1. SI = 0 
Ml == 0 


2. Z ± = Y l - (0.382S + Sl)D i , i = 1, N 

Compute E 
El E 


3. Ml > 0? 


4. Z i = Y i - (0.618S + Sl)D if i * 1, N 

Compute E 
E2 = E 


5. Ml = Ml + 1 


6. Ml = 12? 

7. S = S + 1 

8. E2 > El? 

9. El = E2 

SI = SI + . 382S 
S = . 618S 


10. E2 » El 

S = . 6 1 8S 


Fig. TR-6- 3 . Flow Diagram - Golden Section Search 
For Optimal Step Size 
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Fig. TR-6-4 


Golden Section Search Value Computation Sequence 
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HORIZONTM, PROCESSING 


Introduction 

Because of the motion of the satellite relative to the earth, it 
is a basic property of the swept frequency technique that within 
the same ionogram successive h ' (f ) data points correspond to 
soundings taken through the ionosphere at different earth based 
coordinates (see Fig. TR-7-1) . 

It is implicitly assumed in currently applied data processing 
methods for converting h' (f) data to N (h) profiles that the 
electron density profile is constant over the section of ionosphere 
traversed by the satellite during the period of one ionogram. 

This assumption is required not only for the lamination technique 

r 1 3 i 

of Jackson 1 (forward processing) and Madsen* (inverse processing) 
but also for the least squares estimation technique described 
elsewhere in this report. The error in the electron density 
profile resulting from this assumption is not known. 

Horizontal processing is a method of data reduction which does 
not require the assumption of constant N (h) profile over the 
section required by one ionogram. In horizontal processing 
electron density values along the satellite track are estimated 


*See Section III, this report 
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by linear interpolation between successive ionograms . Horizontal 
processing may be applied in conjunction with either lamination 
or least squares methods of data processing. 

Lamination Method 

In the lamination method it is essential to have a data point 
for the exit frequency so that the electron density at satellite 
height may be obtained. Although we would prefer to linearly 
interpolate electron density at constant height, since the 
satellite height cannot change greatly between successive 
ionograms, we may linearly interpolate to obtain the electron 
density for times between successive ionogram exit frequency 
data points . I t is the interpolated electron density value 
that is used in the computation of electron density for each 
successive lamination. In the Jackson method of Forward Processing 
it is not clear how to select h' (f) points so that the lamination 
altitude will be the same for successive ionograms . With the 
Madsen method of Inverse Processing the lamination boundaries 
may be selected to be the same. All interpolations of electron 
density following the first at satellite height may be made at 
constant height. Of course, extrapolation rather than inter- 
polation will be required for the last ionogram in a sequence . 

The interpolation is illustrated in Fig. TR-7- 2 . At time t,, the 
point 1 is obtained from the exit frequency for ionogram 1. A 
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similar point is obtained for time t 21 from ionogram 2. An h 1 (f) 
data point is selected corresponding to a time t^ 2 » Point 1 is 
now extrapolated forward in time to t^ 2 , giving point 1*. The 
electron density at point 1* is used to compute point 2. This 
process is continued until the last useful data point has been 
processed, say point j. Then we will have a complete N(h) profile 
corresponding to time t^j for the first ionogram, t 2 j for the 
second, etc. 


In the lamination method there are two possible approaches for 
the interpolation of electron density at satellite height. 


The first of these is just to find 


N 


( hij + h^j+i -^lj <hl '3 +1 ' hl 3 > ' t2 :) 




, *23 -hj 

"l.j+l " t lj 


N j+l (h l', j+l ft l, j+l } “ N j (h lj ,t lj ) 


whore 


Nj(h,t) = electron density at height h and time t 

£ . th . 

from 3 ionogram 

th 

h^j * height from i point from j ionogram 
t^j = time of i point from j ionogram 
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The second approach is to find the electron density from the 
ionogram having the higher satellite position at the satellite 
height of the adjacent ionogram without interpolation. This 
step is easily justified because of the very small time interval 
involved in this case. To be specific, suppose that 

h i,j« > h ij 

' Now? using the Madsen Inverse Processing technique, we can find 


N j+l (h l ! , j+l ?t 2, j + l } N j+l (h lj ,t 2, j+l } 

where we have selected h. , . . . = h. . . Using this value of N . , . 

1 ,3+1 I 3 3+1 

we can interpolate at constant height as 

N . (h, . , t„ . ) = N . (h, . , t. .) 

3 13 23 3 13 13 




For either of these cases succeeding interpolations may be made 
at constant height as follows; 


N .. (h . , t, . ) - N.. (h . . , t . . ) 


3 i j k j 


3 1.1 13 


+ 




t. . 

_11 


t ■ ... - t . . 

1,3+1 13 


N.., (h. .,t. N. (h. .,t. .) 

3+l v 13' 1,3+1' 3 13 13 


By this means we may compute successive laminations as shown by 
Fig. TR- 7 - 3 . 
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f~ . , t 
2d 2j 


H * f- 4- 

n 3j' t 3j' t 3j 


h 4 j ' f 4 j ' j 



etc. 


D / • « . . f N 
ionogram number 


Fig. TR-7- 3. 


Horizontal Processing Using Lamination Method 
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Difficulty will be encountered, of course, with ionogram sequences 
in which few successive complete ionograms are available. In 
these cases interpolation in conjunction with a least squares 
estimation may offer a solution. 

Least Squares Estimation 

For this method the initial estimate for N(h) is determined ' 

without extrapolation. Successive revisions to this estimate using 
interpolation are then determined iteratively in order to improve 
the estimate. As shown in Fig. TR-7-4, the initial time skewed N (h) 
profiles are determined without regard to change in N(h) profiles 
with position along the track. Using constant h (or nearly so) 
interpolation, a sequence of N(h) profiles corresponding to 
constant time are found. This sequence of constant time N(h) 
profiles is used in the inverse mapping process in the second 
iteration to find an improved estimate of the time skewed N(h) 
profile. All ionograms in a sequence are processed in this 
manner before proceeding to the next iteration. This process of 
estimating time skewed profiles followed by interpolation may be 
performed sequentially to improve the estimate. As a final step 
the time skewed N(h) profiles are interpolated to the desired 
constant time N(h) profile. 

We know that in actuality the electron density profile is a 
function of position along the satellite trace so that we may 
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write the electron density approximation 


N(h, a ,t) 


For the least squares method we need to use the inverse mapping 
jN.Ch.^tj), i = 1 / • • • • ; 3 j ^ h j ( f j ) t j — 1? • • • • f J 
where J = number of h* (f) data points . 


Since t and f are uniquely related within the time period of one 
ionogram we may write that 


f . <— ■> t . 

3 7 3 


For the first iteration we merely use 


N i (h i? a) , i = 1, 



^ h ( f . ) 
3 3 


However , for subsequent iterations we find 

N . ( h .. , a , t . ) projected on constant 

11 3 time plane 

by 1 i near interpolation as described in the preceding section. 
Therefore all subsequent iterations make use of 

jib (h i , « , t j ) , i = 1, . j} — > h j ( f j ) 


For the least squares method this process is shown in the flow 

diagram of Fig. TR-7- 5. 
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Fig. TPv-7-5. Horizontal Processing Using Least Squares 
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Abstract 

A new computational algorithm is proposed for constructing the 
pseudoinverse of a matrix. The method is believed to be more 
efficient for machine computation than previous methods. 
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Introduction 

A problem which arises frequently in the application of estimation 
theory is to find a best estimate of x from 

Ax = Y- (1) 

where A is m x n 

x is n x 1 
y_ is m x 1 

If A is a square matrix of rank n then x is uniquely determined from 

x = A _1 y (2) 

However, most frequently A is a rectangular matrix with m > n. 

For the case in which the rank of A is n, the least squares estimate 

( 2 ) 

can be found from 

A < T v ’1 T / O \ 

X = (A A) Ay (3) 

since r (A) = n > r (A T A) = n > (A T A) is nonsingular 

(r (A) = rank of A) . 

T 

However , for cases in which r(A) < n it follows that r(A A) < n 
and (A T A) Is singular so that eq . (3) does not apply. For these 
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( X *- 4 } 

cases the pseudoinverse may be used. The pseudoinverse is 

( 4 ) 

defined by the following three requirements : 


I A + Ax = x 


Vx e 

W(A)-L = R (A T ) 

(4) 

* 4 * 

II A z = 0_ 


Vz e 

R (A) = H (A T ) 

(5) 

III A + (£ + z) = 

_ + , , + 
A y; + A _z 

V Z. e 
Vz e 

RCA) = N(A T )-L 
R(A) 1 = W(A T ) 

(6) 

where y e W(A T )~S 

z e N (A T ) 

+ T 

' •> y z 

- 0 

(7) 


M(a)-L is the subspace in E n which is the orthogonal 
complement of the null space of A, 


T T 

R (A ) is the restricted range space of A m E n , 

R (A) is the subspace in which is the orthogonal 
complement of the restricted range of A, 

T T 

M (A ) is the null space of A in E . 

m 

Therefore, we can decompose E and E as follows: 

m n 


R ( A ) 

Cl 

E 

m 

N (A 1 ) 

a 

E 

m 

R (A T ) 

a 

E 

n 

W(A) 

a 

E 

n 
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E m = A/ (A T ) 1 @ W (A T ) 

= R (A) © R (A)"^ 

R (A) = W (A 7 ) 1 

I 

R (A) = W (A 1 ) 

1 

E n = W (A) (+• W (A) 

T X 

= R (A ) (±; R (AT) 

T 1 

R (A ) = N (A) 

m 1 

R (A 1 ) = w (A) 

Notation 

A an m x n matrix 

A ^ inverse of A 

A + pseudoinverse of A 

A pseudo-pseudoinverse of A (meets only requirement I) 

T 

A transpose of A 

W (A) null space of A; z c N (A) <--- > Az = 0_ 

R(A) a restricted range of A 

M (A) orthogonal complement of W (A) 

E space of all column vectors of m elements 

m 

!■' space of all column vectors of n elements 

n 

V for every 

c an element of 

such that 

(+ direct sum of subspaces 

c contained in 
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For the case in which A is a nonsingular square matrix the pseudo- 
inverse becomes the inverse. The pseudoinverse has the property 


that 



( 8 ) 


yields an estimate for x which is best in the sense of leas t 
squares ^ ^ . 

The computational methods described by Greville ^ 1 ^ and Deutsch^ 2 ^ 

( 6 ) 

require the selection of linearly independent sets of rows 
and columns in order to form the pseudoinverse. While this 
approach is very good for an understanding of the construction 
of the pseudoinverse and works well for simple contrived examples, 
it is awkward to apply in practice. 

( 3 ) 

Aoki overcomes this problem by developing a construction of 
the pseudoinverse using eigenvalues and normalized eigenvectors 
of A r A and AA T * . While this method is straightforward and well 
suited for machine computation, the determination of eigenvectors 
is a lengthy process and unattractive from the standpoint of 
efficiency . 


♦Matrices are assumed real so that only simple transpose rather 
than conjugate transpose is used. 
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Proposed Method 


Zadeh and Desoer^^ have shown that 

+ T + T 

A^ = (A'A) A 1 (9) 

so that the determination of the pseudoinverse of any rectangular 
matrix can be reduced to finding the pseudoinverse of a square 
symmetric matrix. 


It is well known that any symmetric n order matrix has n linearly 

independent eigenvectors^ 7 ^. Therefore , according to Pease 
T 

A A is seroisimple and can be written in terms of its eigenvalues 
and projectors as 

A - £ A iPl (10) 

i=l 

T 

where X^ = eigenvalues of A A 

T 

P^ = projectors of A A 


T 

We note that if A A is singular then X^ = 0 for some i, and that 
T 

if r(A A) = r there exists r nonzero values of X^. We define the 
sot S as 

S = li : X. f 0} (11) 


Equation (10) can now be rewritten as 


A A = > X .P . 

1 x 

ieS 
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T, , 


Now we will show that the pseudoinverse of A A is given by 


(A T A) + » £ A.^P. 

ieS 


( 12 ) 


We see that 


(A T A) + (A T A) = Eh p i)( Eh' lp i 
'ieS / 'ieS 


(13) 


( 9 ^ 

Applying the fundamental property of projectors , 


P.P . = 6 . .P. 
i 3 11 i 


(14) 


to equation (13) gives 


(A T A) + (A T A) = £ p. 

ieS 


(15) 


which satisfies the defining relations given in equations (4) , (5) , (6) . 
Therefore, the pseudoinverse is given by eqn. (12) . Next we 
shall develop a recursive method for finding the pseudoinverse. 


Zadeh and Dosoer^ ^ have described Fadeeva's method ^ for finding 
the resolvent (si - A" 1 A) ^ as follows: 


T “1 
(si - A A) 


B(s) 

d(s) 


E 

i=l 


(s 



) 


( 16 ) 
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d(s) = s n + d,s n ^ + d 2 s n ^ 


B (s) 


B, 


B. 


B, 


B, 


n- 1 T 


n-2. 


B 


n~l 

0 


B 0 A T A + d x I 

B^A^A + d2 I 

B, , A T A + d.T 
k-1 k 

B „A T A + d , I 
n-2 n-1 

B ,A T A + d I 
n-1 n 


+ . . . . + d , s + d 
n-r n 

(17) 

« + s B „ + B 

(18) 

n-2 n-l 

d, = -tr(A T A) 

(19) 


2 = -1/2 tr(B x A T A) 

3 = -1/3 tr(B 2 A T A) 

, T 


d k = -1/k tr(B k _ r 

d = -1/n tr (B , A X A) 
n n-1 


A A) 
T, 


where tr 


trace 


The inverse of a nonsingular matrix can be obtained from equations 
(16) , (17) , and (18) as 


T 

(A i A) 


B (0 ) 
d (0 ) 


B 


n-1 


n 


( 20 ) 


However, the pseudoinverse can be found more easily from eqn. 
(19) . 


First, suppose that 

r(A T A) = n (21) 

Since from equation (16) it is clear that the roots of d(s) are 
given by the X ^ , it follows that 

d / 0 
n 


( 22 ) 
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Then from eqn. (19 ) we have 


and 


as before, 


B ,A T A + d I = d> 
n-1 n T 


(A T A) 1 


B 


n-1 


n 


Now suppose that 


r (A T A) = n-1 > 0 


Then since one of the A. - 0 it follows that 

i 


d - 0 
n 


d , f 0 
n-1 


Using these results in eqn. (19) gives 


B ,A T A = (h 
n-1 

(B n A T A + d , I ) A T A = (|) 
n-2 n-1 


Th 


Premultiplying the last equation by A A gives 


R 

m D n .O n p rn m 

(A A A A + A A) A A = <|> 


d 


n-1 


For eqn. (28) to be satisfied 


T ®n~2 T T 

A A — A A + A A = <j> 
d 


(23) 


(24) 


(25) 


(26) 

(27) 


(28) 


( 29 ) 
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The validity of the last statement may be shown as follows. 
From the way in which the are formed in eqn (19) it is 
clear that and A A commute. It is easily verified that 
matrices formed from sums and products of commuting symmetric 
matrices are symmetric*; therefore, from eqn. (19), the B^ 
are symmetric. If two real symmetric matrices commute, then 

( i° ) 

there exists a common diagonalizing matrix, T, such that 


T 1 (A T A)T 


T 1 B.T 

l 


0 X 


n 


v. 


• X 
V i 

n 


(30) 


(31) 


Operating on eqn. (28) with T gives the following sequence 

T _1 (A T A A T A + A T A)A T AT = (j> (32) 

d n-l 


{ (T 1 A T AT) (T 1 - - T) (T 1 A T AT) + T 1 A T AT> (T 1 A T AT) = <j> (33) 

n-1 


* 





~ 

- 

p 



i. 

j 


y 

A 1 

0 1 

X 

j 

i n-2 n 

; v 1 o 



0 


A 1 

0 



A i 

0 



{ 

1 ® 





• 







• 


; # 

: | 





• 




« 


0 

i 

i 

\ 1 
n | 

1 * n-2 

0 V n 


| 

! 0 

A 

n 


0 

A 

n 



0 

* A 

n 


i 

i .J 

L 



- 





T T T 

* (AB) = BA - BA = AB 
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(A . v. n “ 2 A . + A.) A . =0, i = 1, 2, 

'■xi X 11 r 


But r(A 1 A) = n-1 > 0 a n-1 > 0 values of ^ 0 so that 


we must have 


A . v. n 2 X . + A . - * 0 , i = 1 , 2 , ...., n 

11 i x 


m n-2 T T 

< > A A - 3 - - A A + A x A - (|> 

d , 
n-1 


Therefore, eqn. (29) must hold. 


The requirement of eqn. (4) can be restated by premultiplying 
both sides by A to give 

AA + Ax = Ax Vx e M (A) ± = R (A T ) C 


Let us denote by A + ' a matrix which meets requirement (4) but not 
necessarily (5) and (6) . Comparison of eqn. (29) with the 
form of eqn. (4) given above gives 


(A T A) ++ 


We could call A ++ the pseudo-pseudoinverse of A since it satisfies 
only requirement I equation (4) but not necessarily requirements 
II and III, eqn. (5) and eqn. (6). 
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Next we consider the case for which 


r (A 1 A) = n-2 > 0 


Then since two of the = 0 it follows that 


d =0 
n 

d n-l * 0 
d n-2 * 0 


Using these results in eqn» (19) gives 


B _A A = <J> 
n-1 T 


B n-2 (ATA)2 " + 


B , = B ,A T A + d _,I from (19) 
n-z n- Z n z 


(B n _ 3 (A T A) + d n _ 2 I) (A T A) 2 = <f> 


Premultiplying both sides by A A gives 


T ®n- IT T T 2 

(A A - A A + A A) (A A) = <|> 


d 


n-2 


(39) 


(40) 


(41) 

(42) 


(43) 


(44 ) 


For eqn, (44) to be satisfied 
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which can be shown by a 
T 

case for r (A A) = n-1 > 


diagonalizing transformation as in the 
0 . 


Therefore , eqn. (45) must hold so that 


<A T A> ++ = Vi 
n-2 


Continuing in this manner we find that for 


r (A T A) = n-j > 0 


we have 


(A T A) ++ 


B 


n-1- j 


D 


n-3 


(46) 


(4 7 ) 


(48) 


It is not difficult to find examples of matrices such that* 

(A T A) ++ ? (A T A) + (49) 

However, we can show that 

A + = (A T A) ++ A T (50 ) 

is a valid extended version of eqn. (9). First; consider 
requirement I of eqn. (4) in the following form 

I AA + Ax = Ax Vx e N (A)~L 


*A simple 
(A* A) ++ = 


example 
1 

0 1 / 


suggested by C 

do) 


W. Barnes xs A 


1 0 
0 0 


for which 


and A 
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Inserting our proposed solution, given above, yields 


T ++ T 

At (A A) A ] Ax = Ax 

(51) 

.. T, . _ T, , ++_ T,. .It 

A A ( A A ) A Ax ■*“ A Ax 

(52) 


++ 

which is satisfied by the way in whxch A was defined. 


Next we consider 


II A z = 0 


Vz e W(A T ) 


+ 


Inserting our proposed expression for A gives 


T ++ T 
(A X A) A z = 0 


Vz - ' A z = 0 


so that requirement II is satisfied, 


Finally we consider 

III A + (y + z) = A + y_ + A + z Vy_ e W(A T ) 1 

Vz e N(A T ) 

Substituting our proposed solution gives 

(A T A) ++ A T (£ + z) = (A T A) ++ A T y + (A T A) ++ A T z 

so that our proposed solution results in the proper decomposition 
of the space as given by the direct sum 

E n - M (aV + W(A T ) 
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Therefore we have shown that 


4* 


v ++.T 

{A A) A 


T 4*+ 

where (A A) meets only requirement I 


eqn. (4) ) 


Combining this last result with eqn . (48) gives 


+ ® 
a = - a" 

d 

n~D 




(53) 


T 

However , since we are not directly interested in the rank of A A 


we can find A + as 


t 


B 


k - x - a t 


(54) 


where k is the largest value of i for which d. jt 0 


3 

We see from eqn. (19) that approximately n scalar multipli- 
cations are required to evaluate Eh and d . , i = 1» 2, n. 

Therefore, the number of computations rquired to compute the 
pseudoinverse by this method is the same order of magnitude as 
required to invert a nonsingular matrix by the Gauss- Jordan 
method. 


The steps in the computation of the pseudoinverse of A may be 
summarized as follows: 
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T 

1 . Form A A 

2. Compute Bj , , i = 1, 2, . . . . , n from eqn. (19) 

3 . Find k as the largest i for which f 0 

4 . Compute the pseudoinverse of A from eqn . (54 ) 
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Example Problem 

The algorithm described by this paper for finding the pseudo- 
inverse has been programmed in BASIC as shown by the flow 
diagram of Fig. TR-8-1 and associated program listing. This 
program has been tested using the example problem given by 
Greville 1 ^ with results as shown in Table TR-8-1 . These 
results are in agreement with those obtained by Greville . 
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PSEUDO 2-20-70 

t 

LIST 10-490 

SO REM PSEUDO INVERSE PROGRAM 
20 REM AX= V 

30 REM P I S PSEUDO INVERSE 0F A 

40 REM X=PY IS LEAST SQUARES ESTIMATE 

50 REM A HAS M ROWS, N COLUMNS 

•SO REM ENTER "DATA CM),CN>" LINE N0 . 1000 

70 REM ENTER MATRIX ELEMENTS STARTING WITH LINE N0. 1001 

80 REM NUMBER LINES CONSECUTIVELY 
90 REM BEGIN EACH LINE WITH "DATA" 

SOO REM SEPARATE ENTRIES BY COMMAS 

110 REM ENTER ELEMENTS IN FOLLOWING SEQUENCES A BY ROWS, Y 

120 DIM AC10,iO>,XC10>,YClO>,PC10,IO>,BC110,10>,DC10>»CC10,10) 

130 REM B AND D ARE MATRICES IN FADEEVA’S METHOD 

140 REM C-AC TRANSPOSE) *A 

150 READ M, N 

160 FOR 1=1 TO M 

170 FOR J= 1 TO N 

180 READ AC I, J) 

190 NEXT J 

200 NEXT I 

210 F OR 1=1 TO M 

220 READ YCI) 

230 NEXT I 

240 REM FORM C=A C TRANSPOSE ) *A 

250 FOR 1=1 TO N 

260 FOR J=1 TO N 

270 LET CC I, J)=0.0 

280 FOR K= 1 TO M 

290 LET C< I, J)=C< I, J)+A(K, I )*ACK, J) 

300 NEXT K 
310 NEXT J 
320 NEXT I 

330 REM FORM BCN*I+J,K), DC I) BY FADEEVA’S METHOD 

340 LET DC 1 > = 0.0 

350 FOR 1=1 TO N 

360 LET DC 1 ) = DC 1 )-CC I , I ) 

3,70 NEXT 1 

380 F OR J= 1 TO N 

390 FOR K= 1 TO N 

400 LET BCN+J,K)=CC J,K) 

410 NEXT K 

420 NEXT J 

430 FOR J= 1 TO N 

440 LET BCN+J, J) =B CN+ J, J) +D C 1 ) 

450 NEXT J 

460 FOR 1=2 TO N 

470 FOR J= 1 TO N 

480 FOR K= 1 TO N 

490 LET BCN*I+J,K)=0.0 


> 
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PSEUDO 2-20-70 PAGE 2 

t 

LIST 500- 1 100 
500 FOR L= 1 T0 N 

510 LET 8CN*I+J,K)=8CN*I*J J »K)-t-BCN*<I- 1 > •§• Jj> L > *C C L* K } 

520 NEXT L 

530 NEXT K 

540 NEXT J 

550 LET DC I)=0.0 

560 1* OR J=1 TO N 

570 LET DC I > = D< I >-BCN*I +J> J) 

580 NEXT J 

590 LET DC I > = DC I ) / I 

600 FOR J= 1 TO N 

610 LET E3<N*I+J, J) = BCN*I+J, J)+DC I ) 

620 NEXT J 
630 NEXT I 

640 REM FIND LARGEST I SUCH THAT DC I) UNEQUAL TO ZERO 
650 FOR 1=1 TO N 

660 IF ABSCDCN+l-I ) ><1 .E-04 THEN 690 
670 LET I 1 =N+ 1 - I 
680 GOTO 700 
690 NEXT I 

700 REM COMPUTE PSEUDO INVERSE 

710 FOR 1=1 TO N 

720 FOR J=1 TO M 

730 LET P(LJ) = 0.0 

740 FOR K= 1 TO N 

750 LET PCI> J)=PCI, J>-BCN*CI1- 1 )+IjK)*ACJ#K)/DCI 1 ) 

760 NEXT K 

770 NEXT J 

780 NEXT I 

790 REM COMPUTE X 

800 FOR 1=1 TO N 

810 LET XC I >=0.0 

820 FOR J= 1 TO M 

830 LET XC I >=XC I >+PC I> J>*YC J) 

840 NEXT J 
850 NEXT I 

860 PRINT "PSEUDO INVERSE IS" 

870 PRINT 

880 FOR 1=1 TO N 

890 PRINT PCIM)iP(h2)»P(b3> 

900 NEXT I 
910 PRINT 

920 PRINT " I" j "X C I > " 

930 FOR 1=1 TO N 
940 PRINT I i X C I > 

950 NEXT I 
960 PRINT 

970 PRINT "RANK OF A= ” : I 1 
1 ! 00 END 


> 
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L 1ST i 000- i 004 

1000 DATA 3*4 

1001 DATA 4,-1 #-3*2 

1002 DATA -2,5,- 1 *-3 

1003 DATA 2, 1 3, -9, -5 

1004 DATA 7,3,20 
> 



RUN 

PSEUDO INVERSE 

9. 502 969 7fc- 02 
-3 .0790872E-02 
-6.7364R02E-02 
5.0640825E-02 

I 

1 

2 

3 

4 

RANK 0F A= 2 


-5s 6580 1 8 1 E- 02 
3 » 3 1 3535 5E-02 
3. 1 884964E-02 
-3.6730228E-02 

X < I ) 

. 90 l 84433 
• 64035636 
- 1 . 15 73929 
6.6 1 1441 IE -02 


2 .031 885E-02 
3. 782432E-02 
-3.907471 1 E- 02 
-R.9090341L-03 


Table TR-8-1 
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Conclusion 


A computational algorithm using Fadeeva ! s method has been described 
for finding the pseudoinverse. The method is computationally 
efficient in that the number of operations required is of the same 
order as for the Gauss- Jordan method of matrix inversion. 
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APPENDIX INTRODUCTION 


The following eight Technical Reports were written by John Kinkel 
in response to an assignment to investigate the proposed methods 
of closed-loop ionospheric data reduction and provide the mathe- 
matical proofs necessary for a solid theoretical foundation. 

TR-1 "Some Proposed Methods for Reduction of Topside Ionograms 
to Electron Density Profiles" 

Methods are described for finding N(h) or N(h,t) from 
topside ionograms using approximation theory. A gradient 
technique is proposed to find the parameters of the 
approximating function. 

TR-2 "A Weighted Least Squares Approximation Method" 

The theory of Mallinckrodt ' s weighted least squares 
approximation method is investigated. 

TR-3 "Inverse Mapping to Specified f" 

A method is developed for computing values of H(f,a) 
at the same frequency as scaled data points in the h'(f) 
plane . 


A-l 



TR-4 "Simplification of Matrix Inversion Problem" 

A method is developed for solving for X in the matrix 
equation AX = B, without computing A ^ . 

TR-5 "Time Skew Problem" 

The time skew problem in Horizontal Processing is 
investigated. 

TR-6 "Optimal Step Size by One Dimensional Search" 

A method is developed for optimizing the step size of a 
hyper-dimensional adjustment vector that is used in the 
matrix method. 

TR-7 "Horizontal Processing" 

Some of the considerations in Horizontal Processing are 
investigated. 

TR-8 "Computation of the Pseudoinverse" 

The development and proof of a new computational algorithm 
for constructing the pseudoinverse of a matrix is given. 
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